diff --git a/README.md b/README.md
index 143e549..42a5fbf 100644
--- a/README.md
+++ b/README.md
@@ -84,6 +84,7 @@ This option will force Ipopt to use least squares method to calculate dual infea
 
     cy_solver = get_solver(solver="cyipopt-watertap")
     cy_solver.options['recalc_y']='yes'
+    cy_solver.options["recalc_y_feas_tol"] = 1e-2
 
 
 B. Use exact derivatives instead of numeric
diff --git a/src/reaktoro_pse/core/pyomo_property_writer/property_functions.py b/src/reaktoro_pse/core/pyomo_property_writer/property_functions.py
index 13b4615..0c1a808 100644
--- a/src/reaktoro_pse/core/pyomo_property_writer/property_functions.py
+++ b/src/reaktoro_pse/core/pyomo_property_writer/property_functions.py
@@ -10,6 +10,7 @@
 # "https://github.com/watertap-org/reaktoro-pse/"
 #################################################################################
 from pyomo.environ import log10, log, exp
+from idaes.core.util.math import smooth_max
 
 
 def build_scaling_tendency_constraint(rkt_output_object):
@@ -24,6 +25,12 @@ def build_scaling_tendency_constraint(rkt_output_object):
     )
 
 
+# log10(user_output_var) == smooth_max(
+#     build_properties[("saturationIndex", rkt_output_object.property_index)].pyomo_var,
+#     -5,
+# )
+
+
 def build_ph_constraint(rkt_output_object):
     user_output_var = rkt_output_object.pyomo_var
     build_properties = rkt_output_object.pyomo_build_options.properties
diff --git a/src/reaktoro_pse/core/reaktoro_block_builder.py b/src/reaktoro_pse/core/reaktoro_block_builder.py
index 689f8b5..4823567 100644
--- a/src/reaktoro_pse/core/reaktoro_block_builder.py
+++ b/src/reaktoro_pse/core/reaktoro_block_builder.py
@@ -54,14 +54,33 @@ def __init__(self, block, reaktoro_solver, build_on_init=True):
         if isinstance(self.solver, ReaktoroSolver) == False:
             raise TypeError("Reaktoro block builder requires a ReaktoroSolver class")
         self.configure_jacobian_scaling()
+        self.reaktoro_initialize_function = None  # used to provide external solve call
+        self.display_reaktoro_state_function = (
+            None  # used to specifying external function to display rkt state
+        )
+        self.build_output_vars()
         if build_on_init:  # option to support legacy implementation
             self.build_reaktoro_block()
 
-    def build_reaktoro_block(self):
+    def build_reaktoro_block(
+        self,
+        gray_box_model=None,
+        reaktoro_initialize_function=None,
+        display_reaktoro_state_function=None,
+    ):
         """build reaktoro model"""
-        external_model = ReaktoroGrayBox()
-        external_model.configure(self.solver)
-        self.block.reaktoro_model = ExternalGreyBoxBlock(external_model=external_model)
+        if gray_box_model is None:
+            external_model = ReaktoroGrayBox()
+            external_model.configure(self.solver)
+            self.block.reaktoro_model = ExternalGreyBoxBlock(
+                external_model=external_model
+            )
+        else:
+            self.block.reaktoro_model = gray_box_model
+        if reaktoro_initialize_function is not None:
+            self.reaktoro_initialize_function = reaktoro_initialize_function
+        if display_reaktoro_state_function is not None:
+            self.display_reaktoro_state_function = display_reaktoro_state_function
         self.build_input_constraints()
         self.build_output_constraints()
 
@@ -157,31 +176,37 @@ def input_constraints(fs, key):
                         ].get_pyomo_with_required_units()
                     )
 
+    def build_output_vars(self):
+        new_output_vars = {}
+
+        for key, obj in self.solver.output_specs.user_outputs.items():
+            # NOTE: We do not set rkt_outputs to reaktoro_model outputs as they
+            # same as user inputs - we want RKt model to update "user provided vars"
+            # rather then pyomo vars in reaktoro model (e.g. reaktor_block.outputs)
+            if obj.get_pyomo_var() is None:
+                new_output_vars[key] = obj
+        if new_output_vars != {}:
+            self.block.outputs = Var(new_output_vars.keys(), initialize=1)
+            for key, obj in new_output_vars.items():
+                obj.set_pyomo_var(self.block.outputs[key])
+        self.new_output_vars = new_output_vars
+
     def build_output_constraints(self):
         """first update rktOuptutObjects for pyomoBuildProperties with reaktoro pyomo variables as
         they will be used in construction of constraints
         The is will also check if user provided an output pyomo var and if not will
         add them to new_output_var dict, which will be used to create new output variables on the block
         """
-        new_output_vars = {}
         for key, obj in self.solver.output_specs.user_outputs.items():
             if PropTypes.pyomo_built_prop == obj.property_type:
                 for (
                     pyoPropKey,
                     pyoPropObj,
                 ) in obj.pyomo_build_options.properties.items():
-                    pyoPropObj.set_pyomo_var(
-                        self.block.reaktoro_model.outputs[pyoPropKey]
-                    )
-            # NOTE: We do not set rkt_outputs to reaktoro_model outputs as they
-            # same as user inputs - we want RKt model to update "user provided vars"
-            # rather then pyomo vars in reaktoro model (e.g. reaktor_block.outputs)
-            if obj.get_pyomo_var() is None:
-                new_output_vars[key] = obj
-        if new_output_vars != {}:
-            self.block.outputs = Var(new_output_vars.keys(), initialize=1)
-            for key, obj in new_output_vars.items():
-                obj.set_pyomo_var(self.block.outputs[key])
+                    if pyoPropObj.get_pyomo_var() is None:
+                        pyoPropObj.set_pyomo_var(
+                            self.block.reaktoro_model.outputs[pyoPropKey]
+                        )
 
         @self.block.Constraint(self.solver.output_specs.user_outputs)
         def output_constraints(fs, prop, prop_index):
@@ -196,10 +221,14 @@ def output_constraints(fs, prop, prop_index):
                     == self.block.reaktoro_model.outputs[(prop, prop_index)]
                 )
 
-    def initialize(self, presolveDuringInitialization=False):
+    def initialize(self, presolve_during_initialization=False):
         self.initialize_input_variables_and_constraints()
-        self.solver.state.equilibrate_state()
-        self.solver.solve_reaktoro_block(presolve=presolveDuringInitialization)
+
+        if self.reaktoro_initialize_function is None:
+            self.solver.state.equilibrate_state()
+            self.solver.solve_reaktoro_block(presolve=presolve_during_initialization)
+        else:
+            self.reaktoro_initialize_function(presolve=presolve_during_initialization)
         self.initialize_output_variables_and_constraints()
         _log.info(f"Initialized rkt block")
 
@@ -219,9 +248,9 @@ def calc_scale(value):
 
             sf = calc_scale(abs(pyo_var.value))
             if sf > 1e16:
-                _log.warning(f"Var {pyo_var} scale >1e16")
+                _log.warning(f"Var {pyo_var} scale {sf}>1e16")
             if sf < 1e-16:
-                _log.warning(f"Var {pyo_var} scale <1e-16")
+                _log.warning(f"Var {pyo_var} scale {sf}<1e-16")
             return sf
 
     def initialize_output_variables_and_constraints(self):
@@ -267,6 +296,7 @@ def initialize_output_variables_and_constraints(self):
 
         # update jacobian scaling
         self.get_jacobian_scaling()
+        self.set_user_jacobian_scaling()
 
     def get_jacobian_scaling(self):
         if self.jacobian_scaling_type == JacScalingTypes.no_scaling:
@@ -284,12 +314,11 @@ def get_jacobian_scaling(self):
             self.solver.jacobian_scaling_values = (
                 np.sum(np.abs(self.solver.jacobian_matrix) ** 2, axis=1) ** 0.5
             )
-        self.set_user_jacobian_scaling()
 
     def set_user_jacobian_scaling(self, user_scaling=None):
+
         if user_scaling is None:
             user_scaling = self.user_scaling
-
         for i, (key, obj) in enumerate(self.solver.output_specs.rkt_outputs.items()):
             if user_scaling.get(key) != None:
                 scale = user_scaling[key]
@@ -329,3 +358,9 @@ def initialize_input_variables_and_constraints(self):
                 sf = self.get_sf(self.block.reaktoro_model.inputs[key])
             iscale.set_scaling_factor(self.block.reaktoro_model.inputs[key], sf)
             iscale.constraint_scaling_transform(self.block.input_constraints[key], sf)
+
+    def display_state(self):
+        if self.display_reaktoro_state_function is None:
+            print(self.solver.state)
+        else:
+            self.display_reaktoro_state_function()
diff --git a/src/reaktoro_pse/core/reaktoro_gray_box.py b/src/reaktoro_pse/core/reaktoro_gray_box.py
index b18b5f7..8600d26 100644
--- a/src/reaktoro_pse/core/reaktoro_gray_box.py
+++ b/src/reaktoro_pse/core/reaktoro_gray_box.py
@@ -34,13 +34,32 @@ class HessTypes:
 class ReaktoroGrayBox(ExternalGreyBoxModel):
     ########################################################################################
     # custom Grey Box functions
-    def configure(self, reaktoro_solver):
+    def configure(
+        self,
+        reaktoro_solver=None,
+        inputs=None,
+        input_dict=None,
+        outputs=None,
+        hessian_type=None,
+    ):
         # assign a Reaktoro state object to instance
         self.reaktoro_solver = reaktoro_solver
-        self.hess_type = reaktoro_solver.hessian_type
-        self.inputs = reaktoro_solver.input_specs.rkt_inputs.rkt_input_list
-        self.input_dict = reaktoro_solver.input_specs.rkt_inputs
-        self.outputs = list(reaktoro_solver.output_specs.rkt_outputs.keys())
+        if hessian_type is None:
+            self.hess_type = reaktoro_solver.hessian_type
+        else:
+            self.hess_type = hessian_type
+        if inputs is None:
+            self.inputs = reaktoro_solver.input_specs.rkt_inputs.rkt_input_list
+        else:
+            self.inputs = inputs
+        if input_dict is None:
+            self.input_dict = reaktoro_solver.input_specs.rkt_inputs
+        else:
+            self.input_dict = input_dict
+        if outputs is None:
+            self.outputs = list(reaktoro_solver.output_specs.rkt_outputs.keys())
+        else:
+            self.outputs = outputs
         self._outputs_dual_multipliers = np.ones(len(self.outputs))
         self._hess = np.zeros((len(self.inputs), len(self.inputs)))
         self.header_saved = False
@@ -96,7 +115,6 @@ def get_last_output(self, new_params):
                 self.reaktoro_solver.solve_reaktoro_block(params=new_params)
             )
             self.step += 1
-
         self.old_params = copy.deepcopy(new_params)
 
     def evaluate_jacobian_outputs(self):
@@ -112,7 +130,7 @@ def set_output_constraint_multipliers(self, _outputs_dual_multipliers):
         np.copyto(self._outputs_dual_multipliers, _outputs_dual_multipliers)
 
     def get_output_constraint_scaling_factors(self):
-        return self.reaktoro_solver.jacobian_scaling_values
+        return self.reaktoro_solver.get_jacobian_scaling()
 
     def hessian_gauss_newton_version(self, sparse_jac, threshold=7):
 
diff --git a/src/reaktoro_pse/core/reaktoro_inputs.py b/src/reaktoro_pse/core/reaktoro_inputs.py
index 4f0d601..4dd9958 100644
--- a/src/reaktoro_pse/core/reaktoro_inputs.py
+++ b/src/reaktoro_pse/core/reaktoro_inputs.py
@@ -14,6 +14,7 @@
 from reaktoro_pse.core.reaktoro_state import ReaktoroState
 
 import idaes.logger as idaeslog
+import copy
 
 _log = idaeslog.getLogger(__name__)
 
@@ -22,6 +23,44 @@
 """ class to setup input constraints, and specs for reaktoro solver class"""
 
 
+class ReaktoroInputExport:
+    def __init__(self):
+        self.ignore_elements_for_constraints = []
+        self.fixed_solvent_specie = {}
+        self.fixed_solvent_speciation = {}
+        self.rkt_chemical_inputs = None
+        self.assert_charge_neutrality = None
+        self.neutrality_ion = None
+        self.dissolve_species_in_rkt = None
+        self.exact_speciation = None
+
+    def copy_chem_inputs(self, chem_inputs):
+        self.rkt_chemical_inputs = RktInputs()
+        for key, obj in chem_inputs.items():
+            obj.update_values(True)
+            self.rkt_chemical_inputs[key] = None
+            self.rkt_chemical_inputs[key].time_unit = obj.time_unit
+            self.rkt_chemical_inputs[key].main_unit = obj.main_unit
+            self.rkt_chemical_inputs[key].conversion_unit = obj.conversion_unit
+            self.rkt_chemical_inputs[key].conversion_value = obj.conversion_value
+            self.rkt_chemical_inputs[key].required_unit = obj.required_unit
+            self.rkt_chemical_inputs[key].lower_bound = obj.lower_bound
+            self.rkt_chemical_inputs[key].input_type = obj.input_type
+            self.rkt_chemical_inputs[key].value = obj.value
+            self.rkt_chemical_inputs[key].converted_value = obj.converted_value
+        self.rkt_chemical_inputs.registered_phases = chem_inputs.registered_phases
+        self.rkt_chemical_inputs.all_species = chem_inputs.all_species
+        self.rkt_chemical_inputs.species_list = chem_inputs.species_list
+        self.rkt_chemical_inputs.convert_to_rkt_species = (
+            chem_inputs.convert_to_rkt_species
+        )
+        self.rkt_chemical_inputs.composition_is_elements = (
+            chem_inputs.composition_is_elements
+        )
+        self.rkt_chemical_inputs.conversion_method = chem_inputs.conversion_method
+        self.rkt_chemical_inputs.rkt_input_list = chem_inputs.rkt_input_list
+
+
 class ReaktoroInputSpec:
     def __init__(self, reaktor_state):
         # initialize parameters needed to build reaktor solver
@@ -91,7 +130,6 @@ def register_fixed_solvent_specie(self, phase, specie):
         system speciation, so if we want to specify pH, we need to allow system to find eq. H/O and fix
         H2O"""
         self.fixed_solvent_specie[phase] = specie
-
         self.fixed_solvent_speciation[phase] = {}
 
     def register_free_elements(self, elements):
@@ -118,14 +156,16 @@ def configure_specs(
         """
         self.dissolve_species_in_rkt = dissolve_species_in_rkt
         self.exact_speciation = exact_speciation
+
+    def build_input_specs(self):
+        """function to build all the input specs"""
         self.breakdown_species_to_elements()
         self.equilibrium_specs = rkt.EquilibriumSpecs(self.state.state.system())
         self.add_specs(
             self.equilibrium_specs,
             self.assert_charge_neutrality,
-            dissolve_species_in_rkt,
+            self.dissolve_species_in_rkt,
         )
-
         # get input name order!
         for idx, spec in enumerate(self.equilibrium_specs.namesInputs()):
             if spec == "T":
@@ -436,3 +476,29 @@ def write_empty_constraints(self, spec_object):
         for specie in self.empty_constraints:
             spec_object.openTo(specie)
             self.write_empty_con(spec_object, specie)
+
+    def export_config(self):
+        export_object = ReaktoroInputExport()
+        export_object.copy_chem_inputs(self.rkt_chemical_inputs)
+        export_object.ignore_elements_for_constraints = (
+            self.ignore_elements_for_constraints
+        )
+        export_object.fixed_solvent_specie = self.fixed_solvent_specie
+        export_object.fixed_solvent_speciation = self.fixed_solvent_speciation
+        export_object.assert_charge_neutrality = self.assert_charge_neutrality
+        export_object.neutrality_ion = self.neutrality_ion
+        export_object.dissolve_species_in_rkt = self.dissolve_species_in_rkt
+        export_object.exact_speciation = self.exact_speciation
+        return export_object
+
+    def load_from_export_object(self, export_object):
+        self.ignore_elements_for_constraints = (
+            export_object.ignore_elements_for_constraints
+        )
+        self.fixed_solvent_specie = export_object.fixed_solvent_specie
+        self.fixed_solvent_speciation = export_object.fixed_solvent_speciation
+        self.rkt_chemical_inputs = export_object.rkt_chemical_inputs
+        self.assert_charge_neutrality = export_object.assert_charge_neutrality
+        self.neutrality_ion = export_object.neutrality_ion
+        self.dissolve_species_in_rkt = export_object.dissolve_species_in_rkt
+        self.exact_speciation = export_object.exact_speciation
diff --git a/src/reaktoro_pse/core/reaktoro_jacobian.py b/src/reaktoro_pse/core/reaktoro_jacobian.py
index 5e55587..e29896e 100644
--- a/src/reaktoro_pse/core/reaktoro_jacobian.py
+++ b/src/reaktoro_pse/core/reaktoro_jacobian.py
@@ -175,6 +175,13 @@ def display_available(self):
         return available_dict
 
 
+class ReaktoroJacobianExport:
+    def __init__(self):
+        self.der_step_size = None
+        self.jacobian_type = None
+        self.numerical_order = None
+
+
 class ReaktoroJacobianSpec:
     def __init__(self, reaktor_state, reaktor_outputs):
         self.state = reaktor_state
@@ -188,6 +195,20 @@ def __init__(self, reaktor_state, reaktor_outputs):
         self.configure_numerical_jacobian()
         self.check_existing_jacobian_props()
 
+    def export_config(self):
+        export_object = ReaktoroJacobianExport()
+        export_object.der_step_size = self.der_step_size
+        export_object.jacobian_type = self.jacobian_type
+        export_object.numerical_order = self.numerical_order
+        return export_object
+
+    def load_from_export_object(self, export_object):
+        self.configure_numerical_jacobian(
+            jacobian_type=export_object.jacobian_type,
+            order=export_object.numerical_order,
+            step_size=export_object.der_step_size,
+        )
+
     def configure_numerical_jacobian(
         self, jacobian_type="average", order=4, step_size=1e-8
     ):
@@ -200,6 +221,7 @@ def configure_numerical_jacobian(
         """
         self.der_step_size = step_size
         self.jacobian_type = JacType.average
+        self.numerical_order = order
         if jacobian_type == JacType.average:
             self.jacobian_type = JacType.average
             assert order % 2 == 0
diff --git a/src/reaktoro_pse/core/reaktoro_outputs.py b/src/reaktoro_pse/core/reaktoro_outputs.py
index c33851c..32569e9 100644
--- a/src/reaktoro_pse/core/reaktoro_outputs.py
+++ b/src/reaktoro_pse/core/reaktoro_outputs.py
@@ -9,13 +9,17 @@
 # information, respectively. These files are also available online at the URL
 # "https://github.com/watertap-org/reaktoro-pse/"
 #################################################################################
+from matplotlib.pylab import f
 import reaktoro as rkt
 
 import json
 from reaktoro_pse.core.reaktoro_state import ReaktoroState
 import reaktoro_pse.core.pyomo_property_writer.property_functions as propFuncs
 from reaktoro_pse.core.util_classes.rkt_inputs import RktInputTypes
+import copy
+import idaes.logger as idaeslog
 
+_log = idaeslog.getLogger(__name__)
 # disabling warnings
 
 __author__ = "Alexander V. Dudchenko"
@@ -42,13 +46,8 @@ def __init__(
         self.property_type = property_type
         self.property_name = property_name  # properties from which to extract data
         self.property_index = property_index  # index if any
-        if property_type != PropTypes.pyomo_built_prop:
-            # function for getting reaktoro value
-            self.set_get_function(get_function)
-        else:
-            self.set_poyomo_build_option(
-                get_function
-            )  # class that contains information for building pyomo constraints if any
+        self.get_function = None
+        self.set_option_function(property_type, get_function)
         # pyomo var to reference if any - will be built if not user provided
         self.pyomo_var = pyomo_var
         self.value = value  # manually specified value
@@ -67,9 +66,32 @@ def get_value(self, prop_object, update_values=False):
             self.value = value
         return value
 
+    def delete_pyomo_var(self):
+        # self.update_values()
+        del self.pyomo_var
+        self.pyomo_var = None
+
+    def remove_unpicklable_data(self):
+        self.delete_pyomo_var()
+        if self.property_type == PropTypes.pyomo_built_prop:
+            del self.pyomo_build_options
+            # self.get_function = None
+
+    def set_option_function(self, property_type, get_function):
+        if property_type != PropTypes.pyomo_built_prop:
+            # function for getting reaktoro value
+            self.set_get_function(get_function)
+        else:
+            self.set_poyomo_build_option(
+                get_function
+            )  # class that contains information for building pyomo constraints if any
+
     def set_poyomo_build_option(self, func):
         self.pyomo_build_options = func
 
+    def set_pyo_value(self, value):
+        self.pyomo_var.set_value(value)
+
     def set_get_function(self, func):
         self.get_function = func
 
@@ -83,14 +105,19 @@ def get_pyomo_var(self):
         return self.pyomo_var
 
     def set_pyomo_var_value(self, value):
-        self.pyomo_var.value = value
+        self.value = value
+        if self.pyomo_var is not None:
+            self.pyomo_var.set_value(value)
 
-    def get_pyomo_var_value(self, value):
+    def get_pyomo_var_value(self):
         return self.pyomo_var.value
 
     def set_jacobian_value(self, value):
         self.jacobian_value = value
 
+    def get_lb(self):
+        return self.lb
+
 
 class PyomoBuildOptions:
     def __init__(self):
@@ -223,13 +250,47 @@ class PropTypes:
     pyomo_built_prop = "pyomoBuiltProperties"
 
 
+class ReaktoroOutputExport:
+    def __init__(self):
+        self.rkt_outputs = None
+        self.user_outputs = None
+
+    def copy_rkt_outputs(self, outputs):
+        self.rkt_outputs = {}
+        for key, obj in outputs.items():
+            self.rkt_outputs[key] = RktOutput(
+                obj.property_type,
+                obj.property_name,
+                obj.property_index,
+                # get_function=obj.get_function,
+                value=obj.value,
+                jacobian_type=obj.jacobian_type,
+            )
+            self.rkt_outputs[key].remove_unpicklable_data()
+
+    def copy_user_outputs(self, outputs):
+        self.user_outputs = {}  # copy.deepcopy(outputs)
+        for key, obj in outputs.items():
+            self.user_outputs[key] = RktOutput(
+                obj.property_type,
+                obj.property_name,
+                obj.property_index,
+                # get_function=obj.get_function,
+                value=obj.value,
+                jacobian_type=obj.jacobian_type,
+            )
+            self.user_outputs[key].remove_unpicklable_data()
+
+
 class ReaktoroOutputSpec:
     def __init__(self, reaktor_state):
         self.state = reaktor_state
         if isinstance(self.state, ReaktoroState) == False:
             raise TypeError("Reator outputs require rektoroState class")
+
         self.supported_properties = {}
         self.supported_properties[PropTypes.chem_prop] = self.state.state.props()
+
         if RktInputTypes.aqueous_phase in self.state.inputs.registered_phases:
             self.supported_properties[PropTypes.aqueous_prop] = rkt.AqueousProps(
                 self.state.state.props()
@@ -282,6 +343,7 @@ def register_output(
         property_index=None,
         get_all_indexes=False,
         pyomo_var=None,
+        ignore_indexes=None,
     ):
         """register a reaktoro output, couple it to property type.
 
@@ -292,7 +354,7 @@ def register_output(
         pyomo_var -- pyomo var that should be used for the output of this property (optional: will be auto built) (default: None)
         """
         if get_all_indexes:
-            self.get_all_indexes(property_name)
+            self.get_all_indexes(property_name, ignore_indexes)
         else:
             property_type, get_function = self.get_prop_type(
                 property_name, property_index
@@ -313,58 +375,65 @@ def process_output(
         get_function=None,
         pyomo_var=None,
     ):
-        # if property_index is None:
-        #     index = property_name
-        # else:
         index = (property_name, property_index)
-        if property_type != PropTypes.pyomo_built_prop:
-            self.user_outputs[index] = RktOutput(
-                property_type=property_type,
-                property_name=property_name,
-                property_index=property_index,
-                get_function=get_function,
-                pyomo_var=pyomo_var,
-            )
-            if index not in self.rkt_outputs:
-                self.rkt_outputs[index] = self.user_outputs[index]
-        else:
-            self.user_outputs[index] = RktOutput(
-                property_type=property_type,
-                property_name=property_name,
-                property_index=property_index,
-                get_function=get_function,
-                pyomo_var=pyomo_var,
-            )
-            for index, prop in get_function.properties.items():
-                # chcek if prop already exists if it does ont add it outputs
-                # otherwise overwrite it
+        if index not in self.user_outputs:
+            if property_type != PropTypes.pyomo_built_prop:
+                self.user_outputs[index] = RktOutput(
+                    property_type=property_type,
+                    property_name=property_name,
+                    property_index=property_index,
+                    get_function=get_function,
+                    pyomo_var=pyomo_var,
+                )
                 if index not in self.rkt_outputs:
-                    self.rkt_outputs[index] = prop
-                else:
-                    get_function.properties[index] = self.rkt_outputs[index]
+                    self.rkt_outputs[index] = self.user_outputs[index]
+            else:
+                self.user_outputs[index] = RktOutput(
+                    property_type=property_type,
+                    property_name=property_name,
+                    property_index=property_index,
+                    get_function=get_function,
+                    pyomo_var=pyomo_var,
+                )
+                for index, prop in get_function.properties.items():
+                    # chcek if prop already exists if it does nor add it outputs
+                    # otherwise overwrite it
+                    if index not in self.rkt_outputs:
+                        self.rkt_outputs[index] = prop
+                    else:
+                        get_function.properties[index] = self.rkt_outputs[index]
+        else:
+            _log.warning("Output {index}, already added!")
 
     def get_all_indexes(
         self,
         property_name,
+        ignore_indexes,
     ):
         if "species" in property_name:
             for specie in self.species:
-                property_type, get_function = self.get_prop_type(property_name, specie)
-                self.process_output(
-                    property_type=property_type,
-                    property_name=property_name,
-                    property_index=specie,
-                    get_function=get_function,
-                )
+                if ignore_indexes is None or specie not in str(ignore_indexes):
+                    property_type, get_function = self.get_prop_type(
+                        property_name, specie
+                    )
+                    self.process_output(
+                        property_type=property_type,
+                        property_name=property_name,
+                        property_index=specie,
+                        get_function=get_function,
+                    )
         elif "elements" in property_name:
             for element in self.elements:
-                property_type, get_function = self.get_prop_type(property_name, element)
-                self.process_output(
-                    property_type=property_type,
-                    property_name=property_name,
-                    property_index=element,
-                    get_function=get_function,
-                )
+                if ignore_indexes is None or element not in str(ignore_indexes):
+                    property_type, get_function = self.get_prop_type(
+                        property_name, element
+                    )
+                    self.process_output(
+                        property_type=property_type,
+                        property_name=property_name,
+                        property_index=element,
+                        get_function=get_function,
+                    )
         else:
             raise NotImplementedError(
                 f"{property_name} is not supported for automatic indexing"
@@ -405,7 +474,9 @@ def get_prop_type(self, property_name, property_index=None):
                     pass
 
         raise NotImplementedError(
-            f"The {property_name}, {property_index} is not supported at the moment"
+            f"""The {property_name}, {property_index} was not found,
+                its either not supported, or requested index is not in present.
+            """
         )
 
     def get_possible_indexes(self):
@@ -422,6 +493,30 @@ def get_possible_indexes(self):
         #     ].saturationSpecies()
         # ]
 
+    def export_config(self):
+        export_object = ReaktoroOutputExport()
+        export_object.copy_rkt_outputs(self.rkt_outputs)
+        export_object.copy_user_outputs(self.user_outputs)
+        return export_object
+
+    def load_from_export_object(self, export_object):
+        self.rkt_outputs = export_object.rkt_outputs
+        self.user_outputs = export_object.user_outputs
+        for key, obj in self.rkt_outputs.items():
+            property_type, get_function = self.get_prop_type(
+                obj.property_name,
+                obj.property_index,
+            )
+            assert property_type == obj.property_type
+            obj.set_option_function(property_type, get_function)
+        for key, obj in self.user_outputs.items():
+            property_type, get_function = self.get_prop_type(
+                obj.property_name,
+                obj.property_index,
+            )
+            assert property_type == obj.property_type
+            obj.set_option_function(property_type, get_function)
+
     #### start of possible call function to extract values from reactoro properties #####
 
     def _get_prop_phase_name_val(self, prop_type, prop_name, prop_index):
diff --git a/src/reaktoro_pse/core/reaktoro_solver.py b/src/reaktoro_pse/core/reaktoro_solver.py
index bb9787a..2272e06 100644
--- a/src/reaktoro_pse/core/reaktoro_solver.py
+++ b/src/reaktoro_pse/core/reaktoro_solver.py
@@ -33,6 +33,19 @@
 # class to setup reaktor solver for reaktoro
 
 
+class ReaktoroSolverExport:
+    def __init__(self):
+        self.epsilon = None
+        self.tolerance = None
+        self.presolve = None
+        self.presolve_tolerance = None
+        self.presolve_epsilon = None
+        self.max_iters = None
+        self.presolve_max_iters = None
+        self.hessian_type = None
+        self.block_name = None
+
+
 class ReaktoroSolver:
     def __init__(
         self,
@@ -42,7 +55,7 @@ def __init__(
         reaktoro_jacobian_specs,
         block_name=None,
     ):
-        self.blockName = block_name
+        self.block_name = block_name
         self.state = reaktoro_state
         if isinstance(self.state, ReaktoroState) == False:
             raise TypeError("Reator jacobian require rektoroState class")
@@ -52,8 +65,8 @@ def __init__(
         self.output_specs = reaktoro_output_specs
         if isinstance(self.output_specs, ReaktoroOutputSpec) == False:
             raise TypeError("Reator outputs require ReaktoroOutputSpec class")
-        self.jacbian_specs = reaktoro_jacobian_specs
-        if isinstance(self.jacbian_specs, ReaktoroJacobianSpec) == False:
+        self.jacobian_specs = reaktoro_jacobian_specs
+        if isinstance(self.jacobian_specs, ReaktoroJacobianSpec) == False:
             raise TypeError("Reator outputs require ReaktoroOutputSpec class")
         self.solver_options = rkt.EquilibriumOptions()
         self.presolve_options = rkt.EquilibriumOptions()
@@ -74,6 +87,35 @@ def __init__(
         self._sequential_fails = 0
         self._max_fails = 30
         self._input_params = {}
+        self.jacobian_scaling_values = None
+
+    def export_config(self):
+        export_object = ReaktoroSolverExport()
+        export_object.epsilon = self.solver_options.epsilon
+        export_object.tolerance = self.solver_options.optima.convergence.tolerance
+        export_object.presolve = self.presolve
+        export_object.presolve_tolerance = (
+            self.presolve_options.optima.convergence.tolerance
+        )
+        export_object.presolve_epsilon = self.presolve_options.epsilon
+        export_object.max_iters = self.solver_options.optima.maxiters
+        export_object.presolve_max_iters = self.presolve_options.optima.maxiters
+        export_object.hessian_type = self.hessian_type
+        export_object.block_name = self.block_name
+        return export_object
+
+    def load_from_export_object(self, export_object):
+        self.block_name = export_object.block_name
+        self.set_solver_options(
+            export_object.epsilon,
+            export_object.tolerance,
+            export_object.presolve,
+            export_object.presolve_tolerance,
+            export_object.presolve_epsilon,
+            export_object.max_iters,
+            export_object.presolve_max_iters,
+            export_object.hessian_type,
+        )
 
     def set_solver_options(
         self,
@@ -130,6 +172,7 @@ def update_specs(self, params):
             input_obj = self.input_specs.rkt_inputs[input_key]
             if params is None:
                 value = input_obj.get_value(update_temp=True, apply_conversion=True)
+
             else:
                 value = params.get(input_key)
                 input_obj.set_temp_value(value)
@@ -156,11 +199,11 @@ def get_outputs(self):
 
     def get_jacobian(self):
         self.tempJacobianMatrix = self.sensitivity.dudw()
-        self.jacbian_specs.update_jacobian_absolute_values()
+        self.jacobian_specs.update_jacobian_absolute_values()
         jac_matrix = []
         for input_key in self.input_specs.rkt_inputs.rkt_input_list:
             input_obj = self.input_specs.rkt_inputs[input_key]
-            jac_row = self.jacbian_specs.get_jacobian(
+            jac_row = self.jacobian_specs.get_jacobian(
                 self.tempJacobianMatrix, input_obj
             )
             jac_matrix.append(jac_row)
@@ -176,20 +219,25 @@ def solve_reaktoro_block(
         # here we solve reaktor model and return the jacobian matrix and solution, as
         # Cell as update relevant reaktoroSpecs
         self.update_specs(params)
-
-        result = self.try_solve(presolve)
-        self.outputs = self.get_outputs()
-        self.jacobian_matrix = self.get_jacobian()
-        if result.succeeded() == False or display:
+        solve_failed = False
+        try:
+            result = self.try_solve(presolve)
+            self.outputs = self.get_outputs()
+            self.jacobian_matrix = self.get_jacobian()
+        except RuntimeError:
+            solve_failed = True
+        if solve_failed or result.succeeded() == False or display:
             _log.warning(
-                f"warning, solve was not successful for {self.blockName}, fail# {self._sequential_fails}"
+                f"warning, solve was not successful for {self.block_name}, fail# {self._sequential_fails}"
             )
             _log.warning("----inputs were -----")
             for key, value in self._input_params.items():
                 _log.warning(f"{key}: {value}")
             self._sequential_fails += 1
             if self._sequential_fails > self._max_fails:
-                assert False
+                raise RuntimeError(
+                    "Number of failed solves exceed maximum allowed number"
+                )
             raise cyipopt.CyIpoptEvaluationError
         else:
             self._sequential_fails = 0
@@ -214,3 +262,6 @@ def try_solve(self, presolve=False):
         )
         self.output_specs.update_supported_props()
         return result
+
+    def get_jacobian_scaling(self):
+        return self.jacobian_scaling_values
diff --git a/src/reaktoro_pse/core/reaktoro_state.py b/src/reaktoro_pse/core/reaktoro_state.py
index a0a0f1f..a740485 100644
--- a/src/reaktoro_pse/core/reaktoro_state.py
+++ b/src/reaktoro_pse/core/reaktoro_state.py
@@ -9,6 +9,8 @@
 # information, respectively. These files are also available online at the URL
 # "https://github.com/watertap-org/reaktoro-pse/"
 #################################################################################
+import sys
+from numpy import isin
 import reaktoro as rkt
 from pyomo.environ import units as pyunits
 
@@ -16,28 +18,296 @@
 from reaktoro_pse.core.util_classes.rkt_inputs import RktInputTypes
 from pyomo.core.base.var import IndexedVar
 import idaes.logger as idaeslog
+import copy
 
 _log = idaeslog.getLogger(__name__)
 
 __author__ = "Alexander V. Dudchenko"
 
 
-# base class for configuring reaktoro states and solver
+class PhaseData:
+    def __init__(self):
+        """creates state for storing phase object information"""
+        self.phase_list = None
+        self.non_speciate_phase_list = None
+        self.activity_model = None
+        self.state_of_matter = None
+        self.phase_function = None
+        self.phase_list_mode = False
+        self.speciate = True
+
+    def update_species_list(
+        self, species, phase_function, list_mode=False, speciate=True
+    ):
+        """updates list of speciesies for given phase, ensuring no duplication
+
+        args:
+            species - specie or list of specieis to register
+            phase_function - reaktoro function to create phase
+            list_mode - defines if we will need to create a list of phases or single phase
+        """
+
+        def _update_list(current_list, input_val):
+            if isinstance(input_val, str):
+                current_list = [input_val]
+
+            elif isinstance(input_val, list):
+                for i in input_val:
+                    if i not in current_list:
+                        current_list.append(i)
+            else:
+                current_list = input_val
+            return current_list
+
+        if speciate:
+            if species != None and self.phase_list == None:
+                self.phase_list = []
+            self.phase_list = _update_list(self.phase_list, species)
+        else:
+            if species != None and self.non_speciate_phase_list == None:
+                self.non_speciate_phase_list = []
+            self.non_speciate_phase_list = _update_list(
+                self.non_speciate_phase_list, species
+            )
+        self.phase_function = phase_function
+        self.phase_list_mode = list_mode
+
+    def set_activity_model(self, activity_model, state_of_matter=None):
+        """sets activity model and it state"""
+        self.activity_model = activity_model
+        self.state_of_matter = state_of_matter
+
+
+class PhaseManager:
+    def __init__(self):
+        self.registered_phases = {}
+        self.exclude_species_list = []
+
+    def register_phases_species(
+        self, phase, species, phase_function, list_mode=False, speciate=True
+    ):
+        """this is used to add up all phases provided as inputs or registed phases by user,
+        for example, user might provide CO2(g) as a gas input wants to track also N2(g) via register_gas_phase,
+        this ensure we updates all the phases with all requested values
+
+        args:
+            phase - phase type supported by reaktoro-pse
+            species - specie or list of specieis to register
+            phase_function - reaktoro function to create phase
+            list_mode - defines if we will need to create a list of phases or single phase
+        """
+        if phase not in self.registered_phases:
+            self.registered_phases[phase] = PhaseData()
 
+        self.registered_phases[phase].update_species_list(
+            species, phase_function, list_mode, speciate
+        )
 
+    def set_activity_model(
+        self, phase, activity_model, default_activity_model, state_of_matter=None
+    ):
+        """
+        set activity model for reaktoro
+
+        args:
+            phase - phase type supported by reaktoro-pse
+            activity_model - activity model (should be string or touple that contains options
+            to be passed into function)
+            default_activity_model - default activity model if one is not provided
+            state_of_matter - defines state of matter for solids, otherwise none
+        """
+        if phase not in self.registered_phases:
+            self.registered_phases[phase] = PhaseData()
+        if activity_model is None:
+            activity_model = default_activity_model
+        self.registered_phases[phase].set_activity_model(
+            activity_model, state_of_matter
+        )
+
+    def get_registered_phases(self, active_database):
+        """
+        creates listof phases with applied activity models for generation of reaktoro state
+        args:
+            activate_database - database to use during creation of phases
+
+        """
+
+        activate_phase_list = []
+        for phase, phase_object in self.registered_phases.items():
+            if (
+                phase_object.phase_list is not None
+                or phase_object.non_speciate_phase_list is not None
+            ):
+                rkt_phase_object = self.create_rkt_phase(
+                    active_database,
+                    phase_object.phase_function,
+                    phase_object.phase_list,
+                    phase_object.non_speciate_phase_list,
+                    phase_object.phase_list_mode,
+                )
+                if isinstance(rkt_phase_object, list):
+                    for rpo in rkt_phase_object:
+                        self.apply_activity_model(
+                            rpo,
+                            active_database,
+                            phase_object.activity_model,
+                            phase_object.state_of_matter,
+                        )
+                        activate_phase_list.append(rpo)
+                else:
+                    self.apply_activity_model(
+                        rkt_phase_object,
+                        active_database,
+                        phase_object.activity_model,
+                        phase_object.state_of_matter,
+                    )
+                    activate_phase_list.append(rkt_phase_object)
+        return activate_phase_list
+
+    def apply_activity_model(
+        self, rkt_phase_object, active_database, activity_model, state_of_matter=None
+    ):
+        """sets activity mode"""
+        if activity_model is None:
+            raise TypeError(f"Activity model for {rkt_phase_object} is not set.")
+        rkt_activity_model_object = self._process_activity(
+            active_database, activity_model, state_of_matter
+        )
+        rkt_phase_object.set(rkt_activity_model_object)
+
+    def create_rkt_phase(
+        self,
+        active_database,
+        phase_function,
+        phases_list,
+        non_speciate_phase_list,
+        phase_list_mode,
+    ):
+        """Function to remove species from speciation command"""
+        if phase_list_mode:
+            rkt_phase_list = []
+            for phase in phases_list:
+                if isinstance(phase, str):
+                    rkt_phase_list.append(phase_function(phase))
+                else:
+                    rkt_phase_list.append(phase)
+            if non_speciate_phase_list is not None:
+                for phase in non_speciate_phase_list:
+                    if isinstance(phase, str):
+                        rkt_phase_list.append(phase_function(phase))
+                    else:
+                        rkt_phase_list.append(phase)
+            return rkt_phase_list
+        elif isinstance(phases_list, list):
+            if isinstance(phases_list, (list, tuple)):
+                phases_list = " ".join(phases_list)
+            try:
+                # try to spetiatiate
+                phases = phase_function(rkt.speciate(phases_list))
+            except TypeError:
+                # if does nto spectiate, try passing string input directly
+                phases = phase_function(phases_list)
+            system = rkt.ChemicalSystem(active_database, phases)
+            spc_list = []
+            for specie in system.species():
+                if specie.name() not in str(self.exclude_species_list):
+                    spc_list.append(specie.name())
+            if non_speciate_phase_list is not None:
+                for phase in non_speciate_phase_list:
+                    if phase not in spc_list:
+                        spc_list.append(phase)
+            return phase_function(" ".join(spc_list))
+        elif isinstance(non_speciate_phase_list, list):
+            return phase_function(" ".join(non_speciate_phase_list))
+        else:
+            return phases_list
+
+    def _process_activity(
+        self,
+        active_database,
+        activity_model,
+        state_of_matter=None,
+    ):
+        """this will process give activity model if user provides a sting it will
+        find it on reaktoro and initialize it either by it self or passing in state of matter argument
+        if user provides initialized reaktoro activity model we use it directly."""
+
+        def get_func(activity_model, state_of_matter):
+            if isinstance(activity_model, str):
+                if state_of_matter is None:
+                    try:
+                        return getattr(rkt, activity_model)()
+                    except TypeError:
+                        # might require database as input (e.g ActivityModelPhreeqc)
+                        return getattr(rkt, activity_model)(active_database)
+                else:
+                    return getattr(rkt, activity_model)(state_of_matter)
+            elif isinstance(activity_model, tuple):
+                if isinstance(activity_model[1], (tuple, list)):
+                    return getattr(rkt, activity_model[0])(*activity_model[1])
+                else:
+                    return getattr(rkt, activity_model[0])(activity_model[1])
+            else:
+                return activity_model
+
+        if isinstance(activity_model, list):
+            if state_of_matter is None:
+                activity_model = rkt.chain(
+                    *[get_func(am, state_of_matter) for am in activity_model]
+                )
+            else:
+                activity_model = rkt.chain(
+                    *[get_func(am, state_of_matter) for am in activity_model]
+                )
+        else:
+            activity_model = get_func(activity_model, state_of_matter)
+        return activity_model
+
+
+class ReaktoroStateExport:
+    def __init__(self):
+        # defines required exports
+        self.exclude_species_list = None
+        # we can direcly pass original phase manager
+        self.phase_manager = None
+        # this will not have any of the pyomo vars
+        self.inputs = RktInputs.RktInputs()
+        self.exclude_species_list = None
+        self.database_type = None
+        self.database_file = None
+
+    def copy_rkt_inputs(self, inputs):
+        self.inputs = RktInputs.RktInputs()
+        for key, obj in inputs.items():
+            obj.update_values(True)
+            self.inputs[key] = None
+            self.inputs[key].time_unit = obj.time_unit
+            self.inputs[key].main_unit = obj.main_unit
+            self.inputs[key].conversion_unit = obj.conversion_unit
+            self.inputs[key].conversion_value = obj.conversion_value
+            self.inputs[key].required_unit = obj.required_unit
+            self.inputs[key].lower_bound = obj.lower_bound
+            self.inputs[key].input_type = obj.input_type
+            self.inputs[key].value = obj.value
+            self.inputs[key].converted_value = obj.converted_value
+            self.inputs.registered_phases = inputs.registered_phases
+        self.inputs.all_species = inputs.all_species
+        self.inputs.species_list = inputs.species_list
+        self.inputs.convert_to_rkt_species = inputs.convert_to_rkt_species
+        self.inputs.composition_is_elements = inputs.composition_is_elements
+        self.inputs.conversion_method = inputs.conversion_method
+        self.inputs.rkt_input_list = inputs.rkt_input_list
+
+
+# base class for configuring reaktoro states and solver
 class ReaktoroState:
     def __init__(self):
         """initialize all parameters need to build reaktor solver"""
         self.inputs = RktInputs.RktInputs()
-
-        self.mineral_phase = []
-        self.solid_phase = []
-        self.liquid_phase = None
-        self.condensed_phase = None
-        self.aqueous_phase = None
-        self.gas_phase = None
-        self.ion_exchange_phase = None
-        self.registered_phases = {}
+        self.phase_manager = PhaseManager()
+        self.exclude_species_list = []
+        self.export_state = ReaktoroStateExport()
+        self._inputs_not_processed = True
 
     def register_system_inputs(
         self,
@@ -137,6 +407,17 @@ def register_inputs(self, composition, composition_index, phase):
                 self.inputs[specie].set_input_type(phase)
         self._inputs_not_processed = True  # flag that inputs ver modified
 
+    def register_species_to_exclude(self, species):
+        """updates list of species to exclude from when creating phases"""
+        if species != None:
+            if isinstance(species, str):
+                self.phase_manager.exclude_species_list.append(species)
+            elif isinstance(species, list):
+                for spc in species:
+                    self.phase_manager.exclude_species_list.append(spc)
+            else:
+                raise TypeError(f"{species} is not supported, must be str or list")
+
     def register_gas_inputs(
         self,
         composition,
@@ -276,243 +557,139 @@ def process_input(self, var, var_index):
         else:
             return var
 
-    def _process_phase(self, phase, default_phase):
-        """
-        if phsae is list of speices convert it to single string
-        other wise pass it in directly into phase object
-        if phase object is not list, tuple or str assume its user initialized phase
-        and use directly"""
-        if isinstance(phase, (list, tuple)):
-            return default_phase(" ".join(phase))
-        elif isinstance(phase, str):
-            return default_phase(phase)
-        else:
-            return phase
-
-    def _process_phases(self, phase, default_phase, list_mode=False):
-        # List mode is used for when we expect to build a list of phases, such as Mineral and Solid Phase objects
-        if phase != [] and phase is not None:
-            if isinstance(phase, (str)):
-                phase = [phase]
-            if list_mode:
-                phases = []
-                for p in phase:
-                    phases.append(self._process_phase(p, default_phase))
-                return phases
-            else:
-                return self._process_phase(phase, default_phase)
-        else:
-            return None
-
-    def track_registered_phases(self, phase, inputs):
-        """this is used to add up all phases provided as inputs or registed phases by user,
-        for example, user might porbice CO2(g) as a gas input wants to track also N2(g) via register_gas_phase,
-        this ensure we updates all the phases with all requested values"""
-        if inputs is None:
-            return inputs
-        if phase not in self.registered_phases:
-            self.registered_phases[phase] = []
-        if isinstance(self.registered_phases[phase], (str, list)):
-            if isinstance(inputs, str):
-                inputs = [inputs]
-            if isinstance(inputs, list):
-                for i in inputs:
-                    if i not in self.registered_phases[phase]:
-                        self.registered_phases[phase].append(i)
-            else:
-                self.registered_phases[phase] = inputs
-        return self.registered_phases[phase]
-
     def register_aqueous_phase(self, aqueous_phases=None):
         """register possible mineral phases"""
         if aqueous_phases != [] and aqueous_phases is not None:
-            aqueous_phases = self.track_registered_phases(
-                RktInputTypes.aqueous_phase, aqueous_phases
+            self.phase_manager.register_phases_species(
+                RktInputTypes.aqueous_phase, aqueous_phases, rkt.AqueousPhase
             )
-            if isinstance(aqueous_phases, (list, str)):
-                self.aqueous_phase = rkt.AqueousPhase(rkt.speciate(aqueous_phases))
-            else:
-                self.aqueous_phase = aqueous_phases
 
     def register_condensed_phase(self, condensed_phase=None):
         """register possible condensed phases"""
         if condensed_phase != [] and condensed_phase is not None:
-            condensed_phase = self.track_registered_phases(
-                RktInputTypes.condensed_phase, condensed_phase
+            self.phase_manager.register_phases_species(
+                RktInputTypes.condensed_phase, condensed_phase, rkt.CondensedPhases
             )
-            if isinstance(condensed_phase, (list, str)):
-                self.condensed_phase = rkt.CondensedPhases(
-                    rkt.speciate(condensed_phase)
-                )
-            else:
-                self.condensed_phase = condensed_phase
 
     def register_liquid_phase(self, liquid_phases=None):
         """register possible mineral phases"""
         if liquid_phases != [] and liquid_phases is not None:
-            liquid_phases = self.track_registered_phases(
-                RktInputTypes.liquid_phases, liquid_phases
+            self.phase_manager.register_phases_species(
+                RktInputTypes.liquid_phase, liquid_phases, rkt.LiquidPhase
             )
-            if isinstance(liquid_phases, (list, str)):
-                self.liquid_phase = rkt.LiquidPhase(rkt.speciate(liquid_phases))
-            else:
-                self.liquid_phase = liquid_phases
 
     def register_mineral_phases(self, mineral_phases=None):
         """register possible mineral phases"""
-        mineral_phases = self.track_registered_phases(
-            RktInputTypes.mineral_phase, mineral_phases
-        )
-        self.mineral_phase = self._process_phases(
-            mineral_phases, rkt.MineralPhase, list_mode=True
+        self.phase_manager.register_phases_species(
+            RktInputTypes.mineral_phase,
+            mineral_phases,
+            rkt.MineralPhase,
+            list_mode=True,
         )
 
     def register_solid_phases(self, solid_phases=None):
         """register possible solid phases"""
-        solid_phases = self.track_registered_phases(
-            RktInputTypes.solid_phase, solid_phases
-        )
-        self.solid_phase = self._process_phases(
-            solid_phases, rkt.SolidPhase, list_mode=True
+        self.phase_manager.register_phases_species(
+            RktInputTypes.solid_phase, solid_phases, rkt.SolidPhase, list_mode=True
         )
 
-    def register_gas_phase(self, gas_phase=None):
+    def register_gas_phase(self, gas_phase=None, speciate=False):
         """register possible gas phases"""
-        gas_phase = self.track_registered_phases(RktInputTypes.gas_phase, gas_phase)
-        self.gas_phase = self._process_phases(
-            gas_phase, rkt.GaseousPhase, list_mode=False
+        self.phase_manager.register_phases_species(
+            RktInputTypes.gas_phase, gas_phase, rkt.GaseousPhase, speciate=speciate
         )
 
-    def register_ion_exchange_phase(self, ion_exchange_phase=None):
+    def register_ion_exchange_phase(self, ion_exchange_phase=None, speciate=False):
         """register possible ion exchange phases"""
-        ion_exchange_phase = self.track_registered_phases(
-            RktInputTypes.ion_exchange_phase, ion_exchange_phase
-        )
-        self.ion_exchange_phase = self._process_phases(
-            ion_exchange_phase, rkt.IonExchangePhase, list_mode=False
+        self.phase_manager.register_phases_species(
+            RktInputTypes.ion_exchange_phase,
+            ion_exchange_phase,
+            rkt.IonExchangePhase,
+            speciate=speciate,
         )
 
     def set_database(self, dbtype="PhreeqcDatabase", database="pitzer.dat"):
         """set data base of reaktoro"""
         """ assume that if database is string we need to find and init in reaktoro
         other wise assume we received activated reaktoro db"""
-        if isinstance(dbtype, str):
-            self.database = getattr(rkt, dbtype)(database)
+        self.database_type = dbtype
+        self.database_file = database
+
+    def load_database(self):
+        if isinstance(self.database_type, str):
+            self.database = getattr(rkt, self.database_type)(self.database_file)
         else:
-            self.database = dbtype
+            self.database = self.database_type
         self.database_species = [specie.name() for specie in self.database.species()]
         self.database_elements = [
             element.symbol() for element in self.database.elements()
         ]
 
-    def _process_activity(
-        self, activity_model, state_of_matter=None, default_activity_model=None
-    ):
-        """this will process give activity model if user provides a strng it will
-        find it on reaktoro and initialize it either by it self or passing in state of matter argument
-        if user provides intialized reaktoro activity model we use it directly."""
-
-        def get_func(activity_model, state_of_matter):
-            if isinstance(activity_model, str):
-                if state_of_matter is None:
-                    return getattr(rkt, activity_model)()
-                else:
-                    return getattr(rkt, activity_model)(state_of_matter)
-            else:
-                return activity_model
-
-        if activity_model is None:
-            activity_model = default_activity_model
-        if isinstance(activity_model, list):
-            if state_of_matter is None:
-                activity_model = rkt.chain(
-                    *[get_func(am, state_of_matter) for am in activity_model]
-                )
-            else:
-                activity_model = rkt.chain(
-                    *[get_func(am, state_of_matter) for am in activity_model]
-                )
-        else:
-            activity_model = get_func(activity_model, state_of_matter)
-        return activity_model
-
     def set_aqueous_phase_activity_model(
         self,
         activity_model=None,
     ):
         """set activity model of aqueous phases in reaktoro"""
-        self.process_registered_inputs()
-        if self.aqueous_phase is not None:
-            activity_model = self._process_activity(
-                activity_model, default_activity_model="ActivityModelIdealAqueous"
-            )
-            self.aqueous_phase.set(activity_model)
+        self.phase_manager.set_activity_model(
+            RktInputTypes.aqueous_phase,
+            activity_model,
+            default_activity_model="ActivityModelIdealAqueous",
+        )
 
     def set_liquid_phase_activity_model(
         self,
         activity_model=None,
     ):
         """set activity model of liquid phases in reaktoro"""
-        self.process_registered_inputs()
-        if self.liquid_phase is not None:
-            activity_model = self._process_activity(
-                activity_model, default_activity_model="ActivityModelIdealSolution"
-            )
-            self.liquid_phase.set(activity_model)
+
+        self.phase_manager.set_activity_model(
+            RktInputTypes.liquid_phase,
+            activity_model,
+            default_activity_model="ActivityModelIdealSolution",
+        )
 
     def set_gas_phase_activity_model(self, activity_model=None):
         """set activity model of gas phases in reaktoro"""
-        self.process_registered_inputs()
-        activity_model = self._process_activity(
-            activity_model, default_activity_model="ActivityModelIdealGas"
+        self.phase_manager.set_activity_model(
+            RktInputTypes.gas_phase,
+            activity_model,
+            default_activity_model="ActivityModelIdealGas",
         )
-        if self.gas_phase is not None:
-            self.gas_phase.set(activity_model)
 
     def set_condensed_phase_activity_model(self, activity_model=None):
         """set activity model of ccondensed phases in reaktoro"""
-        self.process_registered_inputs()
-        activity_model = self._process_activity(
+        self.phase_manager.set_activity_model(
+            RktInputTypes.condensed_phase,
             activity_model,
             default_activity_model="ActivityModelIdealSolution",
             state_of_matter=rkt.StateOfMatter.Liquid,
         )
-        if self.condensed_phase is not None:
-            self.condensed_phase.set(activity_model)
 
     def set_mineral_phase_activity_model(self, activity_model=None):
         """set activity model of mineral phases in reaktoro"""
-        self.process_registered_inputs()
-        activity_model = self._process_activity(
+        self.phase_manager.set_activity_model(
+            RktInputTypes.mineral_phase,
             activity_model,
             state_of_matter=rkt.StateOfMatter.Solid,
             default_activity_model="ActivityModelIdealSolution",
         )
-        if self.mineral_phase is not None:
-            for phase in self.mineral_phase:
-                phase.set(activity_model)
 
     def set_solid_phase_activity_model(self, activity_model=None):
         """set activity model of solid phases in reaktoro"""
-        self.process_registered_inputs()
-        activity_model = self._process_activity(
+        self.phase_manager.set_activity_model(
+            RktInputTypes.solid_phase,
             activity_model,
             state_of_matter=rkt.StateOfMatter.Solid,
             default_activity_model="ActivityModelIdealSolution",
         )
-        if self.solid_phase is not None:
-            for phase in self.solid_phase:
-                phase.set(activity_model)
 
     def set_ion_exchange_phase_activity_model(self, activity_model=None):
         """set activity model of mineral phases in reaktoro"""
-        self.process_registered_inputs()
-        activity_model = self._process_activity(
+        self.phase_manager.set_activity_model(
+            RktInputTypes.ion_exchange_phase,
             activity_model,
             default_activity_model="ActivityModelIonExchange",
         )
-        if self.ion_exchange_phase is not None:
-            self.ion_exchange_phase.set(activity_model)
 
     def process_registered_inputs(self):
         """this will process all inputs, this function must be called before setting activity models!"""
@@ -537,17 +714,14 @@ def register_all_input_phases(self):
             self.register_aqueous_phase(
                 self.inputs.species_list[RktInputTypes.aqueous_phase]
             )
-
         if len(self.inputs.species_list[RktInputTypes.liquid_phase]) > 0:
             self.register_liquid_phase(
                 self.inputs.species_list[RktInputTypes.liquid_phase]
             )
-
         if len(self.inputs.species_list[RktInputTypes.condensed_phase]) > 0:
             self.register_condensed_phase(
                 self.inputs.species_list[RktInputTypes.condensed_phase]
             )
-
         # register gas phases
         if len(self.inputs.species_list[RktInputTypes.mineral_phase]) > 0:
             self.register_mineral_phases(
@@ -568,20 +742,9 @@ def register_all_input_phases(self):
 
     def build_state(self):
         # this will build reaktor states
-        phases = []
-        if self.aqueous_phase is not None:
-            phases.append(self.aqueous_phase)
-        if self.liquid_phase is not None:
-            phases.append(self.liquid_phase)
-        if self.condensed_phase is not None:
-            phases.append(self.condensed_phase)
-        if self.mineral_phase is not None:
-            for m in self.mineral_phase:
-                phases.append(m)
-        if self.gas_phase is not None:
-            phases.append(self.gas_phase)
-        if self.ion_exchange_phase is not None:
-            phases.append(self.ion_exchange_phase)
+        self.load_database()
+        self.process_registered_inputs()
+        phases = self.phase_manager.get_registered_phases(self.database)
         self.system = rkt.ChemicalSystem(
             self.database,
             *phases,
@@ -631,3 +794,19 @@ def equilibrate_state(self):
         self.set_rkt_state()
         rkt.equilibrate(self.state)
         _log.info("Equilibrated successfully")
+
+    def export_config(self):
+        export_object = ReaktoroStateExport()
+        export_object.copy_rkt_inputs(self.inputs)
+        export_object.exclude_species_list = self.exclude_species_list
+        export_object.phase_manager = self.phase_manager
+        export_object.database_file = self.database_file
+        export_object.database_type = self.database_type
+        return export_object
+
+    def load_from_export_object(self, export_object):
+        self.inputs = export_object.inputs
+        self.phase_manager = export_object.phase_manager
+        self.database_file = export_object.database_file
+        self.database_type = export_object.database_type
+        self.exclude_species_list = export_object.exclude_species_list
diff --git a/src/reaktoro_pse/core/tests/test_reaktoro_builder.py b/src/reaktoro_pse/core/tests/test_reaktoro_builder.py
index 8784004..c9c4d42 100644
--- a/src/reaktoro_pse/core/tests/test_reaktoro_builder.py
+++ b/src/reaktoro_pse/core/tests/test_reaktoro_builder.py
@@ -50,6 +50,7 @@ def build_with_dissolve_in_rkt(build_rkt_state_with_species):
     rkt_inputs.register_chemistry_modifier("CaO", m.lime)
 
     rkt_inputs.configure_specs(dissolve_species_in_rkt=True)
+    rkt_inputs.build_input_specs()
     rkt_outputs = ReaktoroOutputSpec(rkt_state)
     rkt_outputs.register_output("saturationIndex", "Calcite")
     rkt_outputs.register_output("scalingTendency", "Calcite")
@@ -74,6 +75,7 @@ def build_with_dissolve_in_pyomo(build_rkt_state_with_species):
     m.lime.fix()
     rkt_inputs.register_chemistry_modifier("CaO", m.lime)
     rkt_inputs.configure_specs(dissolve_species_in_rkt=False)
+    rkt_inputs.build_input_specs()
     rkt_outputs = ReaktoroOutputSpec(rkt_state)
 
     rkt_outputs.register_output("speciesAmount", get_all_indexes=True)
@@ -97,6 +99,7 @@ def build_with_dissolve_in_rkt_mass_basis(build_rkt_state_with_species_mass_basi
     rkt_inputs.register_chemistry_modifier("CaO", m.lime)
 
     rkt_inputs.configure_specs(dissolve_species_in_rkt=True)
+    rkt_inputs.build_input_specs()
     rkt_outputs = ReaktoroOutputSpec(rkt_state)
     rkt_outputs.register_output("saturationIndex", "Calcite")
     rkt_outputs.register_output("scalingTendency", "Calcite")
@@ -121,6 +124,7 @@ def build_with_dissolve_in_pyomo_mass_basis(build_rkt_state_with_species_mass_ba
     m.lime.fix()
     rkt_inputs.register_chemistry_modifier("CaO", m.lime)
     rkt_inputs.configure_specs(dissolve_species_in_rkt=False)
+    rkt_inputs.build_input_specs()
     rkt_outputs = ReaktoroOutputSpec(rkt_state)
 
     rkt_outputs.register_output("speciesAmount", get_all_indexes=True)
@@ -137,6 +141,8 @@ def test_build_with_rkt_dissolution(build_with_dissolve_in_rkt):
     m.rkt_block = Block()
     builder = ReaktoroBlockBuilder(m.rkt_block, rkt_solver)
     builder.initialize()
+    m.display()
+    m.rkt_block.reaktoro_model.display()
     # will have as many DOFs as outputs due to pyomo not
     # knowing tha graybox exists.
     assert len(m.rkt_block.reaktoro_model.inputs) == len(
@@ -168,6 +174,7 @@ def test_build_with_pyomo_dissolution(build_with_dissolve_in_pyomo):
     builder.initialize()
     # will have as many DOFs as outputs due to pyomo not
     # knowing tha graybox exists.
+    print(rkt_solver.output_specs.rkt_outputs)
     assert degrees_of_freedom(m) == len(rkt_solver.output_specs.rkt_outputs)
     cy_solver = get_solver(solver="cyipopt-watertap")
     cy_solver.options["max_iter"] = 20
diff --git a/src/reaktoro_pse/core/tests/test_reaktoro_inputs.py b/src/reaktoro_pse/core/tests/test_reaktoro_inputs.py
index 6c4af88..9986cf7 100644
--- a/src/reaktoro_pse/core/tests/test_reaktoro_inputs.py
+++ b/src/reaktoro_pse/core/tests/test_reaktoro_inputs.py
@@ -17,6 +17,7 @@
     build_rkt_state_with_species_no_ph,
 )
 
+import pickle
 
 __author__ = "Alexander V. Dudchenko (SLAC)"
 
@@ -29,6 +30,7 @@ def test_with_rkt_sum(build_rkt_state_with_species):
     rkt_state.equilibrate_state()
     rkt_input = RktInputspec.ReaktoroInputSpec(rkt_state)
     rkt_input.configure_specs(dissolve_species_in_rkt=True)
+    rkt_input.build_input_specs()
     input_names = rkt_input.equilibrium_specs.namesControlVariables()
     print(input_names)
     input_constraints = rkt_input.equilibrium_specs.namesConstraints()
@@ -129,6 +131,7 @@ def test_with_rkt_sum_no_ph(build_rkt_state_with_species_no_ph):
     # rkt_input.register_chemistry_modifier("CaO", lime)
     rkt_input.register_charge_neutrality(False)
     rkt_input.configure_specs(dissolve_species_in_rkt=True)
+    rkt_input.build_input_specs()
     input_names = rkt_input.equilibrium_specs.namesControlVariables()
     input_constraints = rkt_input.equilibrium_specs.namesConstraints()
     print("inputtn", input_names)
@@ -201,6 +204,7 @@ def test_with_pyomo_sum(build_rkt_state_with_species):
     rkt_state.equilibrate_state()
     rkt_input = RktInputspec.ReaktoroInputSpec(rkt_state)
     rkt_input.configure_specs(dissolve_species_in_rkt=False)
+    rkt_input.build_input_specs()
     input_names = rkt_input.equilibrium_specs.namesControlVariables()
     input_constraints = rkt_input.equilibrium_specs.namesConstraints()
     print(input_names)
@@ -293,6 +297,7 @@ def test_with_chemical(build_rkt_state_with_species):
 
     rkt_input.register_chemistry_modifier("CaO", lime)
     rkt_input.configure_specs(dissolve_species_in_rkt=True)
+    rkt_input.build_input_specs()
     expected_con_dict = {
         "C": [(1.0, "CO3-2"), (1.0, "CO2")],
         "Na": [(1, "Na+")],
@@ -314,3 +319,46 @@ def test_with_chemical(build_rkt_state_with_species):
         assert eas in rkt_input.active_species
     assert len(rkt_input.active_species) == len(expected_active_species)
     assert len(rkt_input.constraint_dict) == len(expected_con_dict)
+
+
+def test_input_with_pickle_copy(build_rkt_state_with_species):
+    """testing if we can add chemicals"""
+    m, rkt_state = build_rkt_state_with_species
+    rkt_state.build_state()
+    rkt_state.equilibrate_state()
+    rkt_input = RktInputspec.ReaktoroInputSpec(rkt_state)
+
+    lime = Var(initialize=0.01, units=pyunits.mol / pyunits.s)
+    lime.construct()
+
+    rkt_input.register_chemistry_modifier("CaO", lime)
+    rkt_input.configure_specs(dissolve_species_in_rkt=True)
+    rkt_input.build_input_specs()
+    export_object = rkt_input.export_config()
+    pickled_object = pickle.dumps(export_object)
+    unpickled_object = pickle.loads(pickled_object)
+
+    new_rkt_input = RktInputspec.ReaktoroInputSpec(rkt_state)
+    new_rkt_input.load_from_export_object(unpickled_object)
+    new_rkt_input.build_input_specs()
+    expected_con_dict = {
+        "C": [(1.0, "CO3-2"), (1.0, "CO2")],
+        "Na": [(1, "Na+")],
+        "Ca": [(1, "Ca+2"), (1, "CaO")],
+        "Mg": [(1, "Mg+2")],
+    }
+    expected_active_species = ["CaO", "CO3-2", "CO2", "Na+", "Mg+2", "Ca+2"]
+    print(new_rkt_input.constraint_dict)
+    for ion, ecd in expected_con_dict.items():
+        rkt_ecd = new_rkt_input.constraint_dict[ion]
+        # order might change....
+        expected_ecds = len(ecd)
+        counted_ecds = 0
+        for i, (mol, ion) in enumerate(ecd):
+            if mol == rkt_ecd[i][0] and ion == rkt_ecd[i][1]:
+                counted_ecds += 1
+        assert counted_ecds == expected_ecds
+    for eas in expected_active_species:
+        assert eas in new_rkt_input.active_species
+    assert len(new_rkt_input.active_species) == len(expected_active_species)
+    assert len(new_rkt_input.constraint_dict) == len(expected_con_dict)
diff --git a/src/reaktoro_pse/core/tests/test_reaktoro_jacobian.py b/src/reaktoro_pse/core/tests/test_reaktoro_jacobian.py
index 92776e3..859363d 100644
--- a/src/reaktoro_pse/core/tests/test_reaktoro_jacobian.py
+++ b/src/reaktoro_pse/core/tests/test_reaktoro_jacobian.py
@@ -29,6 +29,7 @@
 def build_standard_state(build_rkt_state_with_species):
     m, rkt_state = build_rkt_state_with_species
     rkt_state.register_mineral_phases("Calcite")
+    rkt_state.set_mineral_phase_activity_model()
     rkt_state.build_state()
     rkt_state.equilibrate_state()
     rkt_outputs = ReaktoroOutputSpec(rkt_state)
diff --git a/src/reaktoro_pse/core/tests/test_reaktoro_outputs.py b/src/reaktoro_pse/core/tests/test_reaktoro_outputs.py
index 9c0aba8..5238111 100644
--- a/src/reaktoro_pse/core/tests/test_reaktoro_outputs.py
+++ b/src/reaktoro_pse/core/tests/test_reaktoro_outputs.py
@@ -18,7 +18,7 @@
 from reaktoro_pse.core.tests.test_reaktoro_state import (
     build_rkt_state_with_species,
 )
-
+import pickle
 
 __author__ = "Alexander V. Dudchenko (SLAC)"
 
@@ -85,6 +85,50 @@ def test_output_setup(build_standard_state):
     assert ("pH", None) in rkt_outputs.rkt_outputs
 
 
+def test_output_pickle(build_standard_state, build_rkt_state_with_species):
+    """testing setting out outputs"""
+    rkt_outputs = build_standard_state
+    export_object = rkt_outputs.export_config()
+    pickled_object = pickle.dumps(export_object)
+    unpickled_object = pickle.loads(pickled_object)
+    m, rkt_state = build_rkt_state_with_species
+    new_rkt_outputs = ReaktoroOutputSpec(rkt_state)
+    new_rkt_outputs.load_from_export_object(unpickled_object)
+    """ make sure we have all expected props"""
+    expected_properties = ["chemProp", "aqueousProp", "pyomoBuiltProperties"]
+    for ep in expected_properties:
+        assert ep in rkt_outputs.supported_properties
+    assert len(expected_properties) == len(new_rkt_outputs.supported_properties)
+
+    """test registering chem props"""
+    new_rkt_outputs.register_output("speciesAmount", "Na+")
+    assert ("speciesAmount", "Na+") in new_rkt_outputs.rkt_outputs
+    assert ("speciesAmount", "Na+") in new_rkt_outputs.user_outputs
+    new_rkt_outputs.register_output("speciesAmount", get_all_indexes=True)
+    for specie in new_rkt_outputs.species:
+        assert (
+            new_rkt_outputs.rkt_outputs[("speciesAmount", specie)].property_type
+            == PropTypes.chem_prop
+        )
+    assert len(new_rkt_outputs.rkt_outputs.keys()) == len(new_rkt_outputs.species)
+    new_rkt_outputs.register_output("saturationIndex", "Calcite")
+    assert ("saturationIndex", "Calcite") in new_rkt_outputs.rkt_outputs
+    assert (
+        new_rkt_outputs.rkt_outputs["saturationIndex", "Calcite"].property_type
+        == PropTypes.aqueous_prop
+    )
+    with pytest.raises(NotImplementedError) as a:
+        new_rkt_outputs.register_output("NotRealProp")
+
+    value = new_rkt_outputs.evaluate_property(
+        new_rkt_outputs.rkt_outputs[("speciesAmount", "Na+")]
+    )
+    assert pytest.approx(value, 1e-3) == 0.5
+
+    new_rkt_outputs.register_output("pH")
+    assert ("pH", None) in new_rkt_outputs.rkt_outputs
+
+
 def test_pyomo_constraints(build_standard_state):
     rkt_outputs = build_standard_state
     rkt_outputs.register_output("osmoticPressure", "H2O")
diff --git a/src/reaktoro_pse/core/tests/test_reaktoro_solver.py b/src/reaktoro_pse/core/tests/test_reaktoro_solver.py
index 2419ea8..9919714 100644
--- a/src/reaktoro_pse/core/tests/test_reaktoro_solver.py
+++ b/src/reaktoro_pse/core/tests/test_reaktoro_solver.py
@@ -9,7 +9,14 @@
 # information, respectively. These files are also available online at the URL
 # "https://github.com/watertap-org/reaktoro-pse/"
 #################################################################################
+import enum
+from matplotlib.font_manager import json_load
 import pytest
+
+from reaktoro_pse.core import reaktoro_jacobian
+from reaktoro_pse.core.reaktoro_state import (
+    ReaktoroState,
+)
 from reaktoro_pse.core.reaktoro_jacobian import (
     ReaktoroJacobianSpec,
 )
@@ -26,6 +33,7 @@
 from reaktoro_pse.core.tests.test_reaktoro_state import (
     build_rkt_state_with_species,
 )
+import pickle
 
 __author__ = "Alexander V. Dudchenko (SLAC)"
 
@@ -34,10 +42,12 @@
 def build_standard_state(build_rkt_state_with_species):
     m, rkt_state = build_rkt_state_with_species
     rkt_state.register_mineral_phases("Calcite")
+    rkt_state.set_mineral_phase_activity_model()
     rkt_state.build_state()
     rkt_state.equilibrate_state()
     rkt_inputs = ReaktoroInputSpec(rkt_state)
     rkt_inputs.configure_specs(dissolve_species_in_rkt=True)
+    rkt_inputs.build_input_specs()
     rkt_outputs = ReaktoroOutputSpec(rkt_state)
 
     rkt_outputs.register_output("speciesAmount", get_all_indexes=True)
@@ -54,7 +64,171 @@ def test_solver(build_standard_state):
     rkt_outputs = list(rkt_solver.output_specs.rkt_outputs.keys())
     print(rkt_inputs, rkt_outputs)
     jacobian, outputs = rkt_solver.solve_reaktoro_block()
+    print(rkt_solver.input_specs.constraint_dict)
     print(outputs)
+
+    print([list(l) for l in jacobian])
+    expected_jacboian = [
+        [
+            7.65801494995876e-29,
+            -3.2593494210913563e-36,
+            -2.0741686517690312e-07,
+            -1.2444007911747241e-25,
+            -1.2444007911747241e-25,
+            -6.59081971194882e-25,
+            -4.207269823736039e-25,
+            -1.49021692379296e-25,
+            1.8015999999999957e-09,
+        ],
+        [
+            1.0164395367051604e-20,
+            6.866245319043687e-27,
+            -4.336808689942018e-19,
+            0.0,
+            0.0,
+            7.572474548453445e-18,
+            2.021901526413905e-16,
+            3.608224830031759e-16,
+            1.0000000000000002,
+        ],
+        [
+            8.113696485592391e-09,
+            2.259849376471579e-14,
+            1.8969408354645065e-06,
+            0.00018499351646543246,
+            0.00018499351646543246,
+            -5.365768096869442e-09,
+            -2.079524546710911e-07,
+            -0.00018499351646543243,
+            2.1465127438467364e-08,
+        ],
+        [
+            -2.2978260821996717e-05,
+            -3.5929224848216754e-12,
+            -0.0022255415992898475,
+            0.10861973878704173,
+            0.10861973878704173,
+            2.5854592687684853e-05,
+            -7.862397271260802e-05,
+            -0.1086197387870417,
+            1.2228659892086806e-05,
+        ],
+        [
+            -6.97588777427779e-05,
+            1.752980296255229e-11,
+            -0.004359091473550756,
+            -0.4251074388881366,
+            -0.4251074388881366,
+            1.23303128504726e-05,
+            0.00047786612798534213,
+            0.4251074388881366,
+            6.417220061052413e-05,
+        ],
+        [
+            -9.506954286039246e-05,
+            1.2343205845216623e-11,
+            -0.0067549379919742926,
+            -1.332845072372407,
+            -1.332845072372407,
+            1.0000387004805447,
+            1.9993273694083085,
+            1.332845072372407,
+            7.666185878716834e-05,
+        ],
+        [
+            -4.902071817094089e-05,
+            1.9528117516004683e-11,
+            -0.002301278845038802,
+            0.44991531783507444,
+            0.44991531783507444,
+            -1.3049869556825355e-05,
+            -0.0005057528313723928,
+            -0.44991531783507444,
+            5.220447623445768e-05,
+        ],
+        [
+            -2.319471918733151e-06,
+            -1.5705239528364755e-12,
+            -0.000168062223965199,
+            -0.01617237731406427,
+            -0.01617237731406427,
+            5.078063204763274e-07,
+            0.9989279156972547,
+            0.01617237731406427,
+            2.8189443979070547e-07,
+        ],
+        [
+            2.2319875536741277e-06,
+            1.5720094376045971e-12,
+            0.0001658320299424294,
+            0.016172510973281984,
+            0.016172510973281984,
+            -4.6904451234755533e-07,
+            0.0010624508845250356,
+            -0.016172510973281984,
+            -2.824006434580601e-07,
+        ],
+        [
+            8.748436505902249e-08,
+            -1.4854847681209652e-15,
+            2.2301940227696188e-06,
+            -1.3365921772339677e-07,
+            -1.3365921772339677e-07,
+            -3.876180812891046e-08,
+            9.633418220193148e-06,
+            1.3365921772339675e-07,
+            5.062036675536108e-10,
+        ],
+        [
+            -4.2305648105828074e-23,
+            2.1210189671414445e-28,
+            7.711300019601966e-21,
+            -3.607486416211351e-18,
+            -3.607486416211351e-18,
+            1.0,
+            4.975628662739501e-19,
+            3.469447455833694e-18,
+            2.8879777636005514e-18,
+        ],
+        [
+            4.818680399137594e-09,
+            5.521859127979695e-16,
+            1.3833746784683956e-07,
+            1.4407827958482247e-09,
+            1.4407827958482247e-09,
+            -2.402918241875753e-09,
+            -3.603326006163468e-09,
+            -1.4407827958482247e-09,
+            1.2326277977217798e-09,
+        ],
+        [
+            6.975887774277788e-05,
+            -1.7529802962552286e-11,
+            0.004359091473550756,
+            0.4251074388881365,
+            0.4251074388881365,
+            -1.2330312850461427e-05,
+            -0.0004778661279853019,
+            0.5748925611118636,
+            -6.417220061052484e-05,
+        ],
+        [
+            -8.844850965738781e-18,
+            -6.482170044637456e-13,
+            0.0,
+            -2.592868017854982e-05,
+            -0.0001296434008927491,
+            -5.185736035709962e-07,
+            0.0,
+            2.5928680178549815e-05,
+            -8.271806125530277e-25,
+        ],
+        [0.0, 0.0, 1.0000000002666605, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
+    ]
+    for i, jrow in enumerate(expected_jacboian):
+        for k, jval in enumerate(jrow):
+            assert pytest.approx(jval, 1e-3) == jacobian[i][k]
+
     expected_input_names = [
         "temperature",
         "pressure",
@@ -113,3 +287,262 @@ def test_solver(build_standard_state):
     expected_jac_shape = (len(expected_output_keys), len(expected_input_names))
     assert jacobian.shape[0] == expected_jac_shape[0]
     assert jacobian.shape[1] == expected_jac_shape[1]
+
+
+def test_pickled_solver(build_standard_state):
+    old_rkt_solver = build_standard_state
+    export_state = old_rkt_solver.state.export_config()
+    export_inputs = old_rkt_solver.input_specs.export_config()
+    export_outputs = old_rkt_solver.output_specs.export_config()
+    export_jac = old_rkt_solver.jacobian_specs.export_config()
+    export_solver = old_rkt_solver.export_config()
+
+    export_data = [
+        export_state,
+        export_inputs,
+        export_outputs,
+        export_jac,
+        export_solver,
+    ]
+    pickled_epxort = pickle.dumps(export_data)
+    unpickeld_export = pickle.loads(pickled_epxort)
+
+    rkt_state = ReaktoroState()
+    rkt_state.load_from_export_object(unpickeld_export[0])
+    rkt_state.build_state()
+    rkt_state.equilibrate_state()
+    rkt_inputs = ReaktoroInputSpec(rkt_state)
+    rkt_inputs.load_from_export_object(unpickeld_export[1])
+    rkt_inputs.build_input_specs()
+    print(rkt_inputs.constraint_dict)
+    rkt_outputs = ReaktoroOutputSpec(rkt_state)
+    rkt_outputs.load_from_export_object(unpickeld_export[2])
+    rkt_jacobian = ReaktoroJacobianSpec(rkt_state, rkt_outputs)
+    rkt_jacobian.load_from_export_object(unpickeld_export[3])
+    rkt_solver = ReaktoroSolver(rkt_state, rkt_inputs, rkt_outputs, rkt_jacobian)
+    rkt_solver.load_from_export_object(unpickeld_export[4])
+
+    rkt_inputs = rkt_solver.input_specs.rkt_inputs.rkt_input_list
+    rkt_outputs = list(rkt_solver.output_specs.rkt_outputs.keys())
+
+    jacobian, outputs = rkt_solver.solve_reaktoro_block()
+    print(rkt_solver.state.state)
+    print(outputs)
+    expected_input_names = [
+        "temperature",
+        "pressure",
+        "pH",
+        "CO3-2",
+        "CO2",
+        "Na",
+        "Mg",
+        "Ca",
+        "H2O",
+    ]
+    for i, ein in enumerate(expected_input_names):
+        assert ein in rkt_inputs[i]
+    assert len(expected_input_names) == len(rkt_inputs)
+    expected_output_keys = [
+        ("speciesAmount", "H+"),
+        ("speciesAmount", "H2O"),
+        ("speciesAmount", "CO3-2"),
+        ("speciesAmount", "CO2"),
+        ("speciesAmount", "Ca+2"),
+        ("speciesAmount", "Cl-"),
+        ("speciesAmount", "HCO3-"),
+        ("speciesAmount", "Mg+2"),
+        ("speciesAmount", "MgCO3"),
+        ("speciesAmount", "MgOH+"),
+        ("speciesAmount", "Na+"),
+        ("speciesAmount", "OH-"),
+        ("speciesAmount", "Calcite"),
+        ("saturationIndex", "Calcite"),
+        ("pH", None),
+    ]
+    for i, eon in enumerate(expected_output_keys):
+        assert eon == rkt_outputs[i]
+    assert len(expected_output_keys) == len(rkt_outputs)
+    expected_output_values = [
+        9.00799999999998e-08,
+        50.0,
+        1.2347717588732544e-06,
+        0.0007251176324639623,
+        0.0028374543608618453,
+        0.7024523350480891,
+        0.0030030389116422994,
+        0.09989096781756118,
+        0.00010806304499670852,
+        9.69137442114817e-07,
+        0.5,
+        6.007103894733061e-08,
+        0.007162545639138155,
+        -5.1857360357099634e-15,
+        7.0,
+    ]
+    for i, eov in enumerate(expected_output_values):
+        assert pytest.approx(eov, 1e-3) == outputs[i]
+    assert len(expected_output_values) == len(outputs)
+
+    expected_jac_shape = (len(expected_output_keys), len(expected_input_names))
+    assert jacobian.shape[0] == expected_jac_shape[0]
+    assert jacobian.shape[1] == expected_jac_shape[1]
+    expected_jacboian = [
+        [
+            7.65801494995876e-29,
+            -3.2593494210913563e-36,
+            -2.0741686517690312e-07,
+            -1.2444007911747241e-25,
+            -1.2444007911747241e-25,
+            -6.59081971194882e-25,
+            -4.207269823736039e-25,
+            -1.49021692379296e-25,
+            1.8015999999999957e-09,
+        ],
+        [
+            1.0164395367051604e-20,
+            6.866245319043687e-27,
+            -4.336808689942018e-19,
+            0.0,
+            0.0,
+            7.572474548453445e-18,
+            2.021901526413905e-16,
+            3.608224830031759e-16,
+            1.0000000000000002,
+        ],
+        [
+            8.113696485592391e-09,
+            2.259849376471579e-14,
+            1.8969408354645065e-06,
+            0.00018499351646543246,
+            0.00018499351646543246,
+            -5.365768096869442e-09,
+            -2.079524546710911e-07,
+            -0.00018499351646543243,
+            2.1465127438467364e-08,
+        ],
+        [
+            -2.2978260821996717e-05,
+            -3.5929224848216754e-12,
+            -0.0022255415992898475,
+            0.10861973878704173,
+            0.10861973878704173,
+            2.5854592687684853e-05,
+            -7.862397271260802e-05,
+            -0.1086197387870417,
+            1.2228659892086806e-05,
+        ],
+        [
+            -6.97588777427779e-05,
+            1.752980296255229e-11,
+            -0.004359091473550756,
+            -0.4251074388881366,
+            -0.4251074388881366,
+            1.23303128504726e-05,
+            0.00047786612798534213,
+            0.4251074388881366,
+            6.417220061052413e-05,
+        ],
+        [
+            -9.506954286039246e-05,
+            1.2343205845216623e-11,
+            -0.0067549379919742926,
+            -1.332845072372407,
+            -1.332845072372407,
+            1.0000387004805447,
+            1.9993273694083085,
+            1.332845072372407,
+            7.666185878716834e-05,
+        ],
+        [
+            -4.902071817094089e-05,
+            1.9528117516004683e-11,
+            -0.002301278845038802,
+            0.44991531783507444,
+            0.44991531783507444,
+            -1.3049869556825355e-05,
+            -0.0005057528313723928,
+            -0.44991531783507444,
+            5.220447623445768e-05,
+        ],
+        [
+            -2.319471918733151e-06,
+            -1.5705239528364755e-12,
+            -0.000168062223965199,
+            -0.01617237731406427,
+            -0.01617237731406427,
+            5.078063204763274e-07,
+            0.9989279156972547,
+            0.01617237731406427,
+            2.8189443979070547e-07,
+        ],
+        [
+            2.2319875536741277e-06,
+            1.5720094376045971e-12,
+            0.0001658320299424294,
+            0.016172510973281984,
+            0.016172510973281984,
+            -4.6904451234755533e-07,
+            0.0010624508845250356,
+            -0.016172510973281984,
+            -2.824006434580601e-07,
+        ],
+        [
+            8.748436505902249e-08,
+            -1.4854847681209652e-15,
+            2.2301940227696188e-06,
+            -1.3365921772339677e-07,
+            -1.3365921772339677e-07,
+            -3.876180812891046e-08,
+            9.633418220193148e-06,
+            1.3365921772339675e-07,
+            5.062036675536108e-10,
+        ],
+        [
+            -4.2305648105828074e-23,
+            2.1210189671414445e-28,
+            7.711300019601966e-21,
+            -3.607486416211351e-18,
+            -3.607486416211351e-18,
+            1.0,
+            4.975628662739501e-19,
+            3.469447455833694e-18,
+            2.8879777636005514e-18,
+        ],
+        [
+            4.818680399137594e-09,
+            5.521859127979695e-16,
+            1.3833746784683956e-07,
+            1.4407827958482247e-09,
+            1.4407827958482247e-09,
+            -2.402918241875753e-09,
+            -3.603326006163468e-09,
+            -1.4407827958482247e-09,
+            1.2326277977217798e-09,
+        ],
+        [
+            6.975887774277788e-05,
+            -1.7529802962552286e-11,
+            0.004359091473550756,
+            0.4251074388881365,
+            0.4251074388881365,
+            -1.2330312850461427e-05,
+            -0.0004778661279853019,
+            0.5748925611118636,
+            -6.417220061052484e-05,
+        ],
+        [
+            -8.844850965738781e-18,
+            -6.482170044637456e-13,
+            0.0,
+            -2.592868017854982e-05,
+            -0.0001296434008927491,
+            -5.185736035709962e-07,
+            0.0,
+            2.5928680178549815e-05,
+            -8.271806125530277e-25,
+        ],
+        [0.0, 0.0, 1.0000000002666605, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
+    ]
+    for i, jrow in enumerate(expected_jacboian):
+        for k, jval in enumerate(jrow):
+            assert pytest.approx(jval, 1e-3) == jacobian[i][k]
diff --git a/src/reaktoro_pse/core/tests/test_reaktoro_state.py b/src/reaktoro_pse/core/tests/test_reaktoro_state.py
index b296848..e1da7ad 100644
--- a/src/reaktoro_pse/core/tests/test_reaktoro_state.py
+++ b/src/reaktoro_pse/core/tests/test_reaktoro_state.py
@@ -14,6 +14,7 @@
 import reaktoro_pse.core.reaktoro_state as rktState
 
 from pyomo.environ import ConcreteModel, Var, units as pyunits
+import pickle
 
 __author__ = "Alexander V. Dudchenko (SLAC)"
 
@@ -286,3 +287,21 @@ def test_state_with_elements(build_rkt_state_with_elements):
     assert pytest.approx(float(rkt_state.state.props().elementAmount("Na")), 1e-3) == 0
     assert pytest.approx(float(rkt_state.state.props().elementAmount("Cl")), 1e-3) == 0
     assert pytest.approx(float(rkt_state.state.props().elementAmount("C")), 1e-3) == 0
+
+
+def test_state_with_pickle_copy(build_rkt_state_with_species):
+    m, rkt_state = build_rkt_state_with_species
+    rkt_state.register_mineral_phases("Calcite")
+    rkt_state.set_mineral_phase_activity_model()
+    export_state = rkt_state.export_config()
+    pickeld_state = pickle.dumps(export_state)
+
+    import_state = pickle.loads(pickeld_state)
+    new_rkt_state = rktState.ReaktoroState()
+    new_rkt_state.load_from_export_object(import_state)
+    new_rkt_state.build_state()
+    new_rkt_state.equilibrate_state()
+    assert (
+        pytest.approx(float(new_rkt_state.state.props().speciesAmount("Calcite")), 1e-5)
+        == 0.008835542753970915
+    )
diff --git a/src/reaktoro_pse/core/util_classes/rkt_inputs.py b/src/reaktoro_pse/core/util_classes/rkt_inputs.py
index 7a5477d..a8edb31 100644
--- a/src/reaktoro_pse/core/util_classes/rkt_inputs.py
+++ b/src/reaktoro_pse/core/util_classes/rkt_inputs.py
@@ -60,6 +60,7 @@ def __init__(self, var_name, pyomo_var=None):
         self.main_unit = None
         self.conversion_unit = None
         self.conversion_value = None
+        self.converted_value = None
         self.required_unit = None
         self.rkt_name = var_name
         self.lower_bound = None
@@ -75,6 +76,11 @@ def __init__(self, var_name, pyomo_var=None):
         else:
             self.pyomo_var = None
 
+    def delete_pyomo_var(self):
+        self.update_values(True)
+        del self.pyomo_var
+        self.pyomo_var = None
+
     def get_input_type(self):
         return self.input_type
 
@@ -84,8 +90,10 @@ def set_input_type(self, input_type):
     def update_values(self, update_temp=False):
         if self.pyomo_var is not None:
             self.value = self.pyomo_var.value
-        else:
-            self.value = None
+            if self.conversion_value is not None:
+                self.converted_value = value(self.get_pyomo_with_required_units())
+            else:
+                self.converted_value = self.value
         if update_temp:
             self.set_temp_value(self.value)
 
@@ -104,7 +112,7 @@ def set_lower_bound(self, value):
     def get_value(self, update_temp=False, apply_conversion=False):
         self.update_values(update_temp)
         if apply_conversion:
-            return value(self.get_pyomo_with_required_units())
+            return self.converted_value
         else:
             return self.value
 
@@ -112,12 +120,6 @@ def get_pyomo_with_required_units(self):
         if self.conversion_value == None:
             return self.pyomo_var
         else:
-            # print(
-            #     pyunits.convert(
-            #         self.pyomo_var / (self.conversion_value * self.conversion_unit),
-            #         to_units=self.required_unit,
-            #     )
-            # )
             return pyunits.convert(
                 self.pyomo_var / (self.conversion_value * self.conversion_unit),
                 to_units=self.required_unit,
@@ -168,7 +170,7 @@ def get_rkt_index(self):
         return self.rkt_index
 
     def set_pyomo_var_value(self, value):
-        self.pyomo_var.value = value
+        self.pyomo_var.set_value(value)
 
     def get_pyomo_var_value(self):
         return self.pyomo_var.value
@@ -255,11 +257,16 @@ def _set_species(self, var_name, var, phase):
             if var_name not in self.species_list[phase]:
                 if self.convert_to_rkt_species[phase]:
                     var_name = self.convert_rkt_species_fun(var_name, phase)
-                    super().__setitem__(var_name, var)
-                self.species_list[phase].append(var_name)
+                    if (
+                        var_name not in super().keys()
+                    ):  # make sure its not already thre - can occur if state is reloaded
+                        super().__setitem__(var_name, var)
+                if (
+                    var_name not in self.species_list[phase]
+                ):  # make sure its not already thre - can occur if state is reloaded
+                    self.species_list[phase].append(var_name)
             if var_name not in self.all_species:
                 self.all_species.append(var_name)
-
             if phase not in self.registered_phases:
                 self.registered_phases.append(phase)
 
@@ -275,12 +282,9 @@ def specie_to_rkt_species(species):
     # TODO: needs to be better automated
     name_dict = {
         "-2": ["SO4", "CO3"],
-        "-": [
-            "Cl",
-            "HCO3",
-        ],
+        "-": ["Cl", "HCO3", "F"],
         "+": ["Na", "K"],
-        "+2": ["Mg", "Ca", "Sr", "Ba"],
+        "+2": ["Mg", "Mn", "Ca", "Sr", "Ba"],
         "": ["H2O", "CO2"],
         "H4SiO4": ["Si", "SiO2"],
     }
diff --git a/src/reaktoro_pse/examples/biogas_combustion.py b/src/reaktoro_pse/examples/biogas_combustion.py
index 6c79d13..a761654 100644
--- a/src/reaktoro_pse/examples/biogas_combustion.py
+++ b/src/reaktoro_pse/examples/biogas_combustion.py
@@ -23,6 +23,7 @@
 import idaes.core.util.scaling as iscale
 import reaktoro as rkt
 
+__author__ = "Alexander V. Dudchenko"
 
 # This examples demonstrates how Reaktoro graybox can be used to estimates combustion perofrming an optimization
 # on following reaktoro example: https://reaktoro.org/applications/biomass-gasification/biomass-gasification.html
@@ -128,9 +129,8 @@ def build_biogas(open_species=False):
         gas_phase={
             "composition": m.air,
             "convert_to_rkt_species": False,
-            "phase_components": rkt.GaseousPhase(
-                rkt.speciate("C H O N")
-            ),  # Expects all possible gasses
+            "phase_components": "C H O N",  # Expects all possible gasses
+            "speciate_phase_component": True,
         },
         system_state={
             "pressure": m.feed_pressure,
diff --git a/src/reaktoro_pse/examples/simple_desalination.py b/src/reaktoro_pse/examples/simple_desalination.py
index 0a1adf9..b50dc81 100644
--- a/src/reaktoro_pse/examples/simple_desalination.py
+++ b/src/reaktoro_pse/examples/simple_desalination.py
@@ -27,6 +27,7 @@
 
 import reaktoro as rkt
 
+__author__ = "Alexander V. Dudchenko"
 
 # This examples demonstrates how Reaktoro graybox can be used to estimates
 # properties in desalination process.
diff --git a/src/reaktoro_pse/examples/simple_ion_exchange.py b/src/reaktoro_pse/examples/simple_ion_exchange.py
index 4d79ae9..d309344 100644
--- a/src/reaktoro_pse/examples/simple_ion_exchange.py
+++ b/src/reaktoro_pse/examples/simple_ion_exchange.py
@@ -25,8 +25,7 @@
 import pyomo.environ as pyo
 import reaktoro as rkt
 
-from reaktoro_pse.reaktoro_block_config import jacobian_options
-
+__author__ = "Alexander V. Dudchenko"
 
 # This examples demonstrates how Reaktoro graybox can be used to estimates removal of specific ion through use of ion exchange material.
 
diff --git a/src/reaktoro_pse/examples/tests/test_examples.py b/src/reaktoro_pse/examples/tests/test_examples.py
index 4548cf9..7f748b1 100644
--- a/src/reaktoro_pse/examples/tests/test_examples.py
+++ b/src/reaktoro_pse/examples/tests/test_examples.py
@@ -20,6 +20,7 @@
 
 def test_desal():
     m, m_open = simple_desalination.main()
+
     assert (
         pytest.approx(m.desal_properties[("scalingTendency", "Gypsum")].value, 1e-3)
         == 0.604051223942643
diff --git a/src/reaktoro_pse/examples/thermal_precipitation.py b/src/reaktoro_pse/examples/thermal_precipitation.py
index 11354b7..986b043 100644
--- a/src/reaktoro_pse/examples/thermal_precipitation.py
+++ b/src/reaktoro_pse/examples/thermal_precipitation.py
@@ -24,6 +24,7 @@
 import reaktoro
 import idaes.core.util.scaling as iscale
 
+__author__ = "Alexander V. Dudchenko"
 
 # This examples demonstrates how reaktoro graybox can be used to enthalpy and water vapor pressure.
 
@@ -64,7 +65,7 @@ def build_simple_precipitation():
     m.feed_composition["H2O"].fix(55)
     m.feed_composition["Mg"].fix(0.01)
     m.feed_composition["Na"].fix(0.025)
-    m.feed_composition["Cl"].fix(0.025)
+    m.feed_composition["Cl"].fix(0.03)
     m.feed_composition["Ca"].fix(0.002)
     m.feed_composition["HCO3"].fix(0.01)
     m.feed_composition["SO4"].fix(0.02)
diff --git a/src/reaktoro_pse/parallel_tools/parallel_manager.py b/src/reaktoro_pse/parallel_tools/parallel_manager.py
new file mode 100644
index 0000000..8527f77
--- /dev/null
+++ b/src/reaktoro_pse/parallel_tools/parallel_manager.py
@@ -0,0 +1,360 @@
+#################################################################################
+# WaterTAP Copyright (c) 2020-2024, The Regents of the University of California,
+# through Lawrence Berkeley National Laboratory, Oak Ridge National Laboratory,
+# National Renewable Energy Laboratory, and National Energy Technology
+# Laboratory (subject to receipt of any required approvals from the U.S. Dept.
+# of Energy). All rights reserved.
+#
+# Please see the files COPYRIGHT.md and LICENSE.md for full copyright and license
+# information, respectively. These files are also available online at the URL
+# "https://github.com/watertap-org/reaktoro-pse/"
+#################################################################################
+
+import array
+from calendar import c
+import multiprocessing as mp
+from multiprocessing import shared_memory
+from multiprocessing import Pipe
+from reaktoro_pse.core.reaktoro_state import (
+    ReaktoroState,
+)
+from reaktoro_pse.core.reaktoro_jacobian import (
+    ReaktoroJacobianSpec,
+)
+from reaktoro_pse.core.reaktoro_outputs import (
+    ReaktoroOutputSpec,
+)
+
+from reaktoro_pse.core.reaktoro_inputs import (
+    ReaktoroInputSpec,
+)
+from reaktoro_pse.core.reaktoro_solver import (
+    ReaktoroSolver,
+)
+import numpy as np
+import time
+
+import cyipopt
+import idaes.logger as idaeslog
+
+_log = idaeslog.getLogger(__name__)
+
+
+class RemoteWorker:
+    def __init__(
+        self,
+        config_data,
+        input_reference,
+        output_reference,
+        jacobian_reference,
+    ):
+        """build remote instance of reaktoro state and solver for execution in
+        its own process"""
+        state_config, input_config, output_config, jacobian_config, solver_config = (
+            config_data
+        )
+        self.state = ReaktoroState()
+        self.state.load_from_export_object(state_config)
+        self.state.build_state()
+        self.inputs = ReaktoroInputSpec(self.state)
+        self.inputs.load_from_export_object(input_config)
+        self.inputs.build_input_specs()
+        self.outputs = ReaktoroOutputSpec(self.state)
+        self.outputs.load_from_export_object(output_config)
+        self.jacobian = ReaktoroJacobianSpec(self.state, self.outputs)
+        self.jacobian.load_from_export_object(jacobian_config)
+        self.solver = ReaktoroSolver(
+            self.state, self.inputs, self.outputs, self.jacobian
+        )
+        self.solver.load_from_export_object(solver_config)
+        self.param_dict = {}
+        self.input_reference = shared_memory.SharedMemory(name=input_reference)
+        self.output_reference = shared_memory.SharedMemory(name=output_reference)
+        self.jacobian_reference = shared_memory.SharedMemory(name=jacobian_reference)
+        self.input_matrix = np.ndarray(
+            (3, len(self.inputs.rkt_inputs.keys())),
+            dtype=np.float64,
+            buffer=self.input_reference.buf,
+        )
+        self.old_matrix = None  # np.zeros(len(self.inputs.rkt_inputs.keys()))
+        self.jacobian_matrix = np.ndarray(
+            (
+                len(self.outputs.rkt_outputs.keys()),
+                len(self.inputs.rkt_inputs.keys()),
+            ),
+            dtype=np.float64,
+            buffer=self.jacobian_reference.buf,
+        )
+
+        self.output_matrix = np.ndarray(
+            len(self.outputs.rkt_outputs.keys()),
+            dtype=np.float64,
+            buffer=self.output_reference.buf,
+        )
+        self.params = {}
+        self.old_params = {}
+
+    def initialize(self, presolve=False):
+        _log.info("Initialized in remote worker")
+
+        self.update_inputs()
+        self.state.equilibrate_state()
+        jacobian, outputs = self.solver.solve_reaktoro_block(presolve=presolve)
+        self.update_output_matrix(outputs, jacobian)
+
+    def solve(self):
+        try:
+            if self.check_solve():
+                self.get_params()
+                jacobian, outputs = self.solver.solve_reaktoro_block(self.params)
+                self.update_output_matrix(outputs, jacobian)
+            return WorkerMessages.success
+        except cyipopt.CyIpoptEvaluationError:
+            return WorkerMessages.CyIpoptEvaluationError
+        except RuntimeError:
+            return WorkerMessages.failed
+
+    def update_output_matrix(self, outputs, jacobian):
+        np.copyto(self.output_matrix, outputs)
+        np.copyto(self.jacobian_matrix, jacobian)
+
+    def get_params(self):
+        for i, key in enumerate(self.inputs.rkt_inputs.keys()):
+            self.params[key] = self.input_matrix[2][i]
+
+    def display_state(self):
+        print(self.state.state)
+
+    def check_solve(self):
+        if self.old_matrix is None:
+            self.old_matrix = self.input_matrix[2].copy()
+            return True
+        else:
+            if any(
+                self.old_matrix[i] != self.input_matrix[2][i]
+                for i in range(self.input_matrix[2].size)
+            ):
+                self.old_matrix = self.input_matrix[2].copy()
+                return True
+            else:
+                return False
+
+    def update_inputs(self):
+        for i, key in enumerate(self.inputs.rkt_inputs.keys()):
+            self.inputs.rkt_inputs[key].value = self.input_matrix[0][i]
+            self.inputs.rkt_inputs[key].converted_value = self.input_matrix[1][i]
+
+    def close_shared_memory(self):
+        # clean up memory on termination
+        self.input_reference.close()
+        self.input_reference.unlink()
+        self.output_reference.close()
+        self.output_reference.unlink()
+        self.jacobian_reference.close()
+        self.jacobian_reference.unlink()
+
+
+class WorkerMessages:
+    initialize = "initialize"
+    update_values = "update_values"
+    solve = "solve"
+    success = "success"
+    failed_solve = "failed_solve"
+    CyIpoptEvaluationError = "CyIpoptEvaluationError"
+    terminate = "terminate"
+    failed = "failed"
+    display_state = "display_state"
+
+
+class LocalWorker:
+    def __init__(self, worker_data):
+        """defines local instance of worker to provide
+        direct access to function execution"""
+        self.worker_data = worker_data
+        self.local_pipe, self.remote_pipe = Pipe()
+        self.get_input_and_output_sizes()
+
+    def get_input_and_output_sizes(self):
+        # for storing raw and converted values
+        # index 0 for values for init
+        # index 1 for converted values for init
+        # index 3 for ipopt solver calls
+        self.input_keys = self.worker_data.inputs.rkt_inputs.keys()
+        input_matrix = np.zeros(
+            (3, len(self.worker_data.inputs.rkt_inputs.keys())), dtype=np.float64
+        )
+        self.input_reference = shared_memory.SharedMemory(
+            create=True, size=input_matrix.nbytes
+        )
+        self.input_matrix = np.ndarray(
+            (3, len(self.worker_data.inputs.rkt_inputs.keys())),
+            dtype=np.float64,
+            buffer=self.input_reference.buf,
+        )
+        # for storing output matrix and jacobian
+        jacobian_matrix = np.zeros(
+            (
+                len(self.worker_data.outputs.rkt_outputs.keys()),
+                len(self.input_keys),
+            ),
+            dtype=np.float64,
+        )
+        output_matrix = np.zeros(len(self.worker_data.outputs.rkt_outputs.keys()))
+        self.jacobian_reference = shared_memory.SharedMemory(
+            create=True, size=jacobian_matrix.nbytes
+        )
+        self.jacobian_matrix = np.ndarray(
+            (
+                len(self.worker_data.outputs.rkt_outputs.keys()),
+                len(self.input_keys),
+            ),
+            dtype=np.float64,
+            buffer=self.jacobian_reference.buf,
+        )
+        self.output_reference = shared_memory.SharedMemory(
+            create=True, size=output_matrix.nbytes
+        )
+        self.output_matrix = np.ndarray(
+            len(self.worker_data.outputs.rkt_outputs.keys()),
+            dtype=np.float64,
+            buffer=self.output_reference.buf,
+        )
+
+    def initialize(self, presolve):
+        for i, key in enumerate(self.input_keys):
+
+            self.input_matrix[0][i] = np.float64(
+                self.worker_data.inputs.rkt_inputs[key].get_value()
+            )
+            self.input_matrix[1][i] = np.float64(
+                self.worker_data.inputs.rkt_inputs[key].get_value(apply_conversion=True)
+            )
+        self.local_pipe.send((WorkerMessages.initialize, presolve))
+
+        result = self.local_pipe.recv()
+        # we want to block here.
+
+        if result == WorkerMessages.success:
+            self.update_outputs()
+            _log.warning("Worker initialized")
+        else:
+            raise RuntimeError("Worker failed to initialize")
+
+    def solve(self, params):
+        self.update_params(params)
+        self.local_pipe.send(WorkerMessages.solve)
+
+    def update_outputs(self):
+        for i, key in enumerate(self.worker_data.outputs.rkt_outputs):
+            self.worker_data.outputs.rkt_outputs[key].value = self.output_matrix[i]
+
+    def update_params(self, params):
+        for i, key in enumerate(self.input_keys):
+            self.input_matrix[2][i] = np.float64(params[key])
+
+    def get_solution(self):
+        if self.local_pipe.poll:
+            result = self.local_pipe.recv()
+            if result == WorkerMessages.success:
+                return self.jacobian_matrix.copy(), self.output_matrix.copy()
+            elif result == WorkerMessages.CyIpoptEvaluationError:
+                raise cyipopt.CyIpoptEvaluationError
+            else:
+                raise RuntimeError(
+                    "The worker failed and did not return a solution terminated"
+                )
+
+    def display_state(self):
+        self.local_pipe.send(WorkerMessages.display_state)
+
+    def terminate(self):
+        self.local_pipe.send(WorkerMessages.terminate)
+        _log.info("Worker terminated")
+
+
+class ReaktoroParallelManager:
+    def __init__(self, time_out):
+        self.registered_workers = {}
+        self.processes = {}
+        self.time_out = time_out
+
+    def register_block(self, block_idx, block_data):
+        self.registered_workers[block_idx] = LocalWorker(block_data)
+
+    def get_solve_and_get_function(self, block_idx):
+        return (
+            self.registered_workers[block_idx].solve,
+            self.registered_workers[block_idx].get_solution,
+        )
+
+    def get_initialize_function(self, block_idx):
+        return self.registered_workers[block_idx].initialize
+
+    def get_display_function(self, block_idx):
+        return self.registered_workers[block_idx].display_state
+
+    def start_workers(self):
+        for idx, local_worker in self.registered_workers.items():
+            process = mp.Process(
+                target=ReaktoroActor,
+                args=(
+                    local_worker.remote_pipe,
+                    local_worker.worker_data.frozen_state,
+                    local_worker.input_reference.name,
+                    local_worker.output_reference.name,
+                    local_worker.jacobian_reference.name,
+                    self.time_out,
+                ),
+            )
+            process.start()
+            _log.info(f"Started parallel worker {idx}")
+            self.processes[idx] = process
+
+    def terminate_workers(self):
+        for idx, local_worker in self.registered_workers.items():
+            local_worker.terminate()
+            self.processes[idx].join()
+
+
+def ReaktoroActor(
+    pipe, reaktoro_block_data, input_matrix, output_matrix, jacobian_matrix, time_out=20
+):
+    reaktoro_worker = RemoteWorker(
+        reaktoro_block_data, input_matrix, output_matrix, jacobian_matrix
+    )
+    dog_watch = time.time()
+    while True:
+        if pipe.poll():
+            msg = pipe.recv()
+            if isinstance(msg, tuple):
+                command = msg[0]
+                option = msg[1]
+            else:
+                command = msg
+
+            if command == WorkerMessages.update_values:
+                reaktoro_worker.update_inputs()
+                result = WorkerMessages.success
+            if command == WorkerMessages.initialize:
+                reaktoro_worker.initialize(presolve=option)
+                result = WorkerMessages.success
+            if command == WorkerMessages.solve:
+                result = reaktoro_worker.solve()
+            if command == WorkerMessages.display_state:
+                result = reaktoro_worker.display_state()
+                result = WorkerMessages.success
+            if command == WorkerMessages.terminate:
+                reaktoro_worker.close_shared_memory()
+                return
+            pipe.send(result)
+            dog_watch = time.time()
+        if abs(time.time() - dog_watch) > time_out:
+            # make sure we kill worker if it does not receive command with in time out
+            # this is to handle scenario where main flowsheet crashes or is rebuilt
+            _log.warning(
+                f"""Worker timed out, shutting down worker.
+                    The time out was set to {time_out}, increase it if necessary in ReaktoroBlockManager options"""
+            )
+            reaktoro_worker.close_shared_memory()
+            return
+        time.sleep(1e-3)  # 1ms sleep time to reduce load when idle
diff --git a/src/reaktoro_pse/parallel_tools/reaktoro_block_manager.py b/src/reaktoro_pse/parallel_tools/reaktoro_block_manager.py
new file mode 100644
index 0000000..084ae02
--- /dev/null
+++ b/src/reaktoro_pse/parallel_tools/reaktoro_block_manager.py
@@ -0,0 +1,332 @@
+#################################################################################
+# WaterTAP Copyright (c) 2020-2024, The Regents of the University of California,
+# through Lawrence Berkeley National Laboratory, Oak Ridge National Laboratory,
+# National Renewable Energy Laboratory, and National Energy Technology
+# Laboratory (subject to receipt of any required approvals from the U.S. Dept.
+# of Energy). All rights reserved.
+#
+# Please see the files COPYRIGHT.md and LICENSE.md for full copyright and license
+# information, respectively. These files are also available online at the URL
+# "https://github.com/watertap-org/reaktoro-pse/"
+#################################################################################
+
+import multiprocessing
+from matplotlib.pyplot import sca
+from pytest import param
+from reaktoro_pse.core.reaktoro_gray_box import (
+    ReaktoroGrayBox,
+)
+from reaktoro_pse.parallel_tools.parallel_manager import ReaktoroParallelManager
+
+from reaktoro_pse.core.reaktoro_gray_box import HessTypes
+
+from pyomo.contrib.pynumero.interfaces.external_grey_box import (
+    ExternalGreyBoxBlock,
+)
+import numpy as np
+from idaes.core.base.process_base import declare_process_block_class, ProcessBlockData
+from pyomo.common.config import ConfigValue, IsInstance, ConfigDict
+
+
+class ReaktoroBlockData:
+    def __init__(self):
+        self.state = None
+        self.inputs = None
+        self.outputs = None
+        self.solver = None
+        self.jacobian = None
+        self.builder = None
+        self.pseudo_gray_box = None
+
+    def get_configs(self):
+        configs = []
+        configs.append(self.state.export_config())
+        configs.append(self.inputs.export_config())
+        configs.append(self.outputs.export_config())
+        configs.append(self.jacobian.export_config())
+        configs.append(self.solver.export_config())
+        return configs
+
+    def freeze_state(self):
+        self.frozen_state = self.get_configs()
+
+
+class AggregateSolverState:
+    def __init__(self, parallel_mode=True, maximum_number_of_parallel_solves=None):
+        self.hessian_type = None
+        self.inputs = []
+        self.input_dict = {}
+        self.input_blk_indexes = []
+        self.registered_blocks = None
+        self.outputs = []
+        self.output_blk_indexes = []
+        self.jacobian_scaling_obj = []
+        self.solver_functions = {}
+        self.get_solution_function = {}
+        self.output_matrix = []
+        self.jacobian_matrix = []
+        self.input_windows = {}
+        self.output_windows = {}
+        self.parallel_mode = parallel_mode
+        if maximum_number_of_parallel_solves is None:
+            maximum_number_of_parallel_solves = multiprocessing.cpu_count() - 1
+        self.maximum_number_of_parallel_solves = maximum_number_of_parallel_solves
+
+    def register_solve_function(self, block_index, solver_function):
+        self.solver_functions[block_index] = solver_function
+
+    def register_get_function(self, block_index, get_function):
+        self.get_solution_function[block_index] = get_function
+
+    def register_input(self, block_index, input_key, input_obj):
+        self.inputs.append((block_index, input_key))
+        self.input_dict[(block_index, input_key)] = input_obj
+        self.input_dict[(block_index, input_key)].original_key = input_key
+        self.input_blk_indexes.append(block_index)
+
+    def register_output(self, block_index, output_key):
+        self.outputs.append((block_index, output_key))
+        self.output_blk_indexes.append(block_index)
+
+        if len(self.output_blk_indexes) > 1:
+            self.jacobian_matrix = np.zeros((len(self.outputs), len(self.inputs)))
+            self.output_matrix = np.zeros(len(self.outputs))
+            self.get_windows(block_index)
+
+    def get_windows(self, block_idx):
+        _, output_unique_sets = np.unique(self.output_blk_indexes, return_inverse=True)
+        self.registered_blocks, input_unique_sets = np.unique(
+            self.input_blk_indexes, return_inverse=True
+        )
+        input_idx = np.arange(len(self.inputs))
+        output_idx = np.arange(len(self.outputs))
+        output_start, output_end = min(
+            output_idx[output_unique_sets == block_idx]
+        ), max(output_idx[output_unique_sets == block_idx])
+        input_start, input_end = min(input_idx[input_unique_sets == block_idx]), max(
+            input_idx[input_unique_sets == block_idx]
+        )
+        self.input_windows[block_idx] = (
+            input_start,
+            input_end + 1,
+        )  # need to offset by 1
+        self.output_windows[block_idx] = (
+            output_start,
+            output_end + 1,
+        )  # need to offset by 1
+
+    def register_scaling_vals(self, scaling_values):
+        self.jacobian_scaling_obj.append(scaling_values)
+
+    def get_jacobian_scaling(self):
+        scaling_array = None
+        for obj in self.jacobian_scaling_obj:
+            if scaling_array is None:
+                scaling_array = obj
+            else:
+                np.hstack(scaling_array, obj)
+        return scaling_array
+
+    def get_params(self, block_idx, params):
+        param_set = {}
+        for (idx, key), item in params.items():
+            if block_idx == idx:
+                param_set[key] = item
+        return param_set
+
+    def update_solution(self, block_idx, output, jacobian):
+        self.output_matrix[
+            self.output_windows[block_idx][0] : self.output_windows[block_idx][1]
+        ] = output
+        self.jacobian_matrix[
+            self.output_windows[block_idx][0] : self.output_windows[block_idx][1],
+            self.input_windows[block_idx][0] : self.input_windows[block_idx][1],
+        ] = jacobian
+
+    def solve_reaktoro_block(self, params):
+        if self.parallel_mode:
+            return self.parallel_solver(params)
+        else:
+            return self.serial_solver(params)
+
+    def parallel_solver(self, params):
+        active_workers = []
+        for blk in self.registered_blocks:
+            param_set = self.get_params(blk, params)
+            self.solver_functions[blk](param_set)
+            active_workers.append(blk)
+            # if we have more than max workers,
+            # collect any results that are ready first
+            while len(active_workers) >= self.maximum_number_of_parallel_solves:
+                for blk in active_workers:
+                    result = self.get_solution_function[active_workers[0]]()
+                    if result is not None:
+                        jacobian, output = result
+                        self.update_solution(blk, output, jacobian)
+                        active_workers.pop(0)  # remove first worker
+                        break
+        # collect any results that are still not collected
+        for blk in active_workers:
+            jacobian, output = self.get_solution_function[blk]()
+            self.update_solution(blk, output, jacobian)
+
+        return (
+            self.jacobian_matrix,
+            self.output_matrix,
+        )
+
+    def serial_solver(self, params):
+        for blk in self.registered_blocks:
+            param_set = self.get_params(blk, params)
+            jacobian, output = self.solver_functions[blk](param_set)
+            self.update_solution(blk, output, jacobian)
+        return (
+            self.jacobian_matrix,
+            self.output_matrix,
+        )
+
+
+class PseudoGrayBox:
+    def __init__(self):
+        self.inputs = {}
+        self.outputs = {}
+
+    def register_input(self, aggregate_inputs, input_keys, block_idx):
+        for input_key in input_keys:
+            self.inputs[input_key] = aggregate_inputs[(block_idx, input_key)]
+
+    def register_output(self, aggregate_outputs, output_keys, block_idx):
+        for output_keys in output_keys:
+            self.outputs[output_keys] = aggregate_outputs[(block_idx, output_keys)]
+
+
+@declare_process_block_class("ReaktoroBlockManager")
+class ReaktoroBlockManagerData(ProcessBlockData):
+    CONFIG = ProcessBlockData.CONFIG()
+    CONFIG.declare(
+        "hessian_type",
+        ConfigValue(
+            default="BFGS",
+            domain=IsInstance((str, HessTypes)),
+            description="Hessian type to use for reaktor gray box",
+            doc="""Hessian type to use, some might provide better stability
+                options (Jt.J, BFGS, BFGS-mod, BFGS-damp, BFGS-ipopt""",
+        ),
+    )
+    CONFIG.declare(
+        "use_parallel_mode",
+        ConfigValue(
+            default=True,
+            domain=bool,
+            description="Enables use of parallel workers",
+            doc="""If true, will parallelize all rekatoro solver calls using multiprocessing""",
+        ),
+    )
+    CONFIG.declare(
+        "worker_timeout",
+        ConfigValue(
+            default=20,
+            domain=int,
+            description="Defines time in seconds for worker time out",
+            doc="""This is time out for parallel workers to time out and shut down if they receive no 
+            commands from main process (e.g. flowsheet)""",
+        ),
+    )
+    CONFIG.declare(
+        "maximum_number_of_parallel_solves",
+        ConfigValue(
+            default=None,
+            domain=int,
+            description="Maximum number of parallel solves",
+            doc="""This will limit how many parallel solves would be run, this will not reduce number of 
+            spawned processes""",
+        ),
+    )
+
+    def build(self):
+        super().build()
+        self.registered_blocks = []
+        self.aggregate_solver_state = AggregateSolverState(
+            parallel_mode=self.config.use_parallel_mode,
+            maximum_number_of_parallel_solves=self.config.maximum_number_of_parallel_solves,
+        )
+        if self.config.use_parallel_mode:
+            self.parallel_manager = ReaktoroParallelManager(
+                self.config.worker_timeout,
+            )
+
+    def register_block(self, state, inputs, outputs, jacobian, solver, builder):
+        blk = ReaktoroBlockData()
+        blk.state = state
+        blk.inputs = inputs
+        blk.outputs = outputs
+        blk.jacobian = jacobian
+        blk.solver = solver
+        blk.builder = builder
+        blk.freeze_state()
+
+        self.registered_blocks.append(blk)
+        return blk
+
+    def aggregate_inputs_and_outputs(self):
+        for block_idx, block in enumerate(self.registered_blocks):
+            if self.config.use_parallel_mode:
+                self.parallel_manager.register_block(block_idx, block)
+            for key, obj in block.inputs.rkt_inputs.items():
+                self.aggregate_solver_state.register_input(block_idx, key, obj)
+            for output in block.outputs.rkt_outputs.keys():
+                self.aggregate_solver_state.register_output(block_idx, output)
+            if self.config.use_parallel_mode:
+                solve_func, get_func = self.parallel_manager.get_solve_and_get_function(
+                    block_idx
+                )
+                self.aggregate_solver_state.register_solve_function(
+                    block_idx, solve_func
+                )
+                self.aggregate_solver_state.register_get_function(block_idx, get_func)
+            else:
+                self.aggregate_solver_state.register_solve_function(
+                    block_idx, block.solver.solve_reaktoro_block
+                )
+        if self.config.use_parallel_mode:
+            self.parallel_manager.start_workers()
+
+    def build_reaktoro_blocks(self):
+        self.aggregate_inputs_and_outputs()
+        external_model = ReaktoroGrayBox()
+        external_model.configure(
+            self.aggregate_solver_state,
+            inputs=self.aggregate_solver_state.inputs,
+            input_dict=self.aggregate_solver_state.input_dict,
+            outputs=self.aggregate_solver_state.outputs,
+            hessian_type=self.config.hessian_type,
+        )
+        self.reaktoro_model = ExternalGreyBoxBlock(external_model=external_model)
+        for block_idx, block in enumerate(self.registered_blocks):
+            pseudo_gray_box_model = PseudoGrayBox()
+            pseudo_gray_box_model.register_input(
+                self.reaktoro_model.inputs,
+                block.inputs.rkt_inputs.keys(),
+                block_idx,
+            )
+            pseudo_gray_box_model.register_output(
+                self.reaktoro_model.outputs,
+                block.outputs.rkt_outputs.keys(),
+                block_idx,
+            )
+            if self.config.use_parallel_mode:
+                init_func = self.parallel_manager.get_initialize_function(block_idx)
+                disp_func = self.parallel_manager.get_display_function(block_idx)
+            else:
+                init_func = self.aggregate_solver_state.solver_functions[block_idx]
+                disp_func = None
+            block.builder.build_reaktoro_block(
+                gray_box_model=pseudo_gray_box_model,
+                reaktoro_initialize_function=init_func,
+                display_reaktoro_state_function=disp_func,
+            )
+            block.pseudo_gray_box = pseudo_gray_box_model
+
+    def terminate_workers(self):
+        if self.config.use_parallel_mode:
+            self.parallel_manager.terminate_workers()
diff --git a/src/reaktoro_pse/parallel_tools/tests/test_manager.py b/src/reaktoro_pse/parallel_tools/tests/test_manager.py
new file mode 100644
index 0000000..1a34170
--- /dev/null
+++ b/src/reaktoro_pse/parallel_tools/tests/test_manager.py
@@ -0,0 +1,171 @@
+#################################################################################
+# WaterTAP Copyright (c) 2020-2024, The Regents of the University of California,
+# through Lawrence Berkeley National Laboratory, Oak Ridge National Laboratory,
+# National Renewable Energy Laboratory, and National Energy Technology
+# Laboratory (subject to receipt of any required approvals from the U.S. Dept.
+# of Energy). All rights reserved.
+#
+# Please see the files COPYRIGHT.md and LICENSE.md for full copyright and license
+# information, respectively. These files are also available online at the URL
+# "https://github.com/watertap-org/reaktoro-pse/"
+#################################################################################
+
+import pytest
+
+from reaktoro_pse.parallel_tools.reaktoro_block_manager import ReaktoroBlockManager
+
+from reaktoro_pse.reaktoro_block import ReaktoroBlock
+from reaktoro_pse.tests.test_reaktoro_block import (
+    build_rkt_state_with_species,
+)
+from pyomo.environ import (
+    ConcreteModel,
+    Var,
+    assert_optimal_termination,
+    units as pyunits,
+)
+from watertap.core.solvers import get_solver
+
+
+def test_blockBuild_with_speciation_block(build_rkt_state_with_species):
+    m = build_rkt_state_with_species
+    m.CaO = Var(["CaO"], initialize=0.001, units=pyunits.mol / pyunits.s)
+    m.CaO.fix()
+    m.reaktoro_manager = ReaktoroBlockManager()
+    m.property_block = ReaktoroBlock(
+        aqueous_phase={
+            "composition": m.composition,
+            "convert_to_rkt_species": True,
+        },
+        system_state={
+            "temperature": m.temp,
+            "pressure": m.pressure,
+            "pH": m.pH,
+        },
+        jacobian_options={
+            "numerical_type": "average",
+            "numerical_order": 2,
+            "numerical_step": 1e-8,
+        },
+        database="PhreeqcDatabase",
+        database_file="pitzer.dat",
+        chemistry_modifier=m.CaO,
+        outputs=m.outputs,
+        build_speciation_block=True,
+        reaktoro_block_manager=m.reaktoro_manager,
+    )
+    m.reaktoro_manager.build_reaktoro_blocks()
+    m.property_block.initialize()
+
+    cy_solver = get_solver(solver="cyipopt-watertap")
+    cy_solver.options["max_iter"] = 20
+    m.pH.unfix()
+    m.outputs[("scalingTendency", "Calcite")].fix(5)
+    result = cy_solver.solve(m, tee=True)
+    assert_optimal_termination(result)
+    m.display()
+    assert pytest.approx(m.outputs[("pH", None)].value, 1e-2) == 6.7496301
+    assert pytest.approx(m.pH.value, 1e-2) == 6.401
+
+    m.property_block.display_jacobian_outputs()
+
+    scaling_result = m.property_block.display_jacobian_scaling()
+    expected_scaling = {
+        "speciation_block": {
+            ("speciesAmount", "H+"): 9.007999999999993e-08,
+            ("speciesAmount", "H2O"): 50.0,
+            ("speciesAmount", "CO3-2"): 3.2175702176273733e-06,
+            ("speciesAmount", "CO2"): 0.00189035577659813,
+            ("speciesAmount", "Ca+2"): 0.01,
+            ("speciesAmount", "Cl-"): 0.7116050981506346,
+            ("speciesAmount", "HCO3-"): 0.007825323588838813,
+            ("speciesAmount", "Mg+2"): 0.09971792990850152,
+            ("speciesAmount", "MgCO3"): 0.0002811030643454316,
+            ("speciesAmount", "MgOH+"): 9.670271530541402e-07,
+            ("speciesAmount", "Na+"): 0.5,
+            ("speciesAmount", "OH-"): 6.004424745615723e-08,
+        },
+        "property_block": {
+            ("saturationIndex", "Calcite"): 1.554873983061197,
+            ("pH", None): 7.520409745594153,
+        },
+    }
+    assert "speciation_block" in scaling_result
+    assert "property_block" in scaling_result
+    new_scaling = {}
+    for key in scaling_result["speciation_block"]:
+        new_scaling[key] = 1
+        assert (
+            pytest.approx(scaling_result["speciation_block"][key], 1e-3)
+            == expected_scaling["speciation_block"][key]
+        )
+    m.property_block.update_jacobian_scaling(new_scaling)
+    scaling_result = m.property_block.display_jacobian_scaling()
+
+    print(m.property_block.reaktoro_model.outputs)
+    assert len(m.property_block.reaktoro_model.outputs) == 2
+
+    assert "speciation_block" in scaling_result
+    for key in scaling_result["speciation_block"]:
+        assert scaling_result["speciation_block"][key] == 1
+    new_scaling = {}
+    for key in scaling_result["property_block"]:
+        new_scaling[key] = 1
+        assert (
+            pytest.approx(scaling_result["property_block"][key], 1e-3)
+            == expected_scaling["property_block"][key]
+        )
+    m.property_block.update_jacobian_scaling(new_scaling)
+    scaling_result = m.property_block.display_jacobian_scaling()
+
+    assert "property_block" in scaling_result
+    for key in scaling_result["property_block"]:
+        assert scaling_result["property_block"][key] == 1
+    m.property_block.display_reaktoro_state()
+    m.reaktoro_manager.terminate_workers()
+
+
+def test_blockBuild_with_wateqf_data_base(build_rkt_state_with_species):
+    m = build_rkt_state_with_species
+    m.CaO = Var(["CaO"], initialize=0.001, units=pyunits.mol / pyunits.s)
+    m.CaO.fix()
+    m.reaktoro_manager = ReaktoroBlockManager()
+    m.property_block = ReaktoroBlock(
+        aqueous_phase={
+            "composition": m.composition,
+            "convert_to_rkt_species": True,
+            "activity_model": "ActivityModelPhreeqc",
+        },
+        system_state={
+            "temperature": m.temp,
+            "pressure": m.pressure,
+            "pH": m.pH,
+        },
+        database="PhreeqcDatabase",
+        database_file="wateq4f.dat",
+        chemistry_modifier=m.CaO,
+        outputs=m.outputs,
+        build_speciation_block=True,
+        reaktoro_block_manager=m.reaktoro_manager,
+        exclude_species_list=[
+            "CH4",
+            "O2",
+            "S2-2",
+            "S4-2",
+            "S3-2",
+            "S5-2",
+            "S6-2",
+        ],
+    )
+    m.reaktoro_manager.build_reaktoro_blocks()
+    m.property_block.initialize()
+    cy_solver = get_solver(solver="cyipopt-watertap")
+    cy_solver.options["max_iter"] = 20
+    m.pH.unfix()
+    m.outputs[("scalingTendency", "Calcite")].fix(5)
+    result = cy_solver.solve(m, tee=True)
+    assert_optimal_termination(result)
+    m.display()
+    assert pytest.approx(m.outputs[("pH", None)].value, 1e-2) == 7.49301431889365
+    assert pytest.approx(m.pH.value, 1e-2) == 6.515501990042
+    m.reaktoro_manager.terminate_workers()
diff --git a/src/reaktoro_pse/reaktoro_block.py b/src/reaktoro_pse/reaktoro_block.py
index 2b8a535..99194e2 100644
--- a/src/reaktoro_pse/reaktoro_block.py
+++ b/src/reaktoro_pse/reaktoro_block.py
@@ -31,6 +31,7 @@
 from reaktoro_pse.core.reaktoro_solver import ReaktoroSolver
 from reaktoro_pse.core.reaktoro_block_builder import ReaktoroBlockBuilder
 
+from reaktoro_pse.parallel_tools.reaktoro_block_manager import ReaktoroBlockManager
 
 from reaktoro_pse.reaktoro_block_config.jacobian_options import JacobianOptions
 from reaktoro_pse.reaktoro_block_config.reaktoro_solver_options import (
@@ -61,7 +62,10 @@ class ReaktoroBlockData(ProcessBlockData):
             include_pH=True, aqueous_phase=False, include_solvent_species=True
         ),
     )
-    CONFIG.declare(RktInputTypes.gas_phase, PhaseInput().get_dict())
+    CONFIG.declare(
+        RktInputTypes.gas_phase,
+        PhaseInput().get_dict(include_speciate_phase_component=True),
+    )
     CONFIG.declare(RktInputTypes.condensed_phase, PhaseInput().get_dict())
     CONFIG.declare(RktInputTypes.mineral_phase, PhaseInput().get_dict())
     CONFIG.declare(RktInputTypes.solid_phase, PhaseInput().get_dict())
@@ -79,6 +83,29 @@ class ReaktoroBlockData(ProcessBlockData):
              property state is modified through addition of chemicals or formation of phases that modify final state""",
         ),
     )
+    CONFIG.declare(
+        "exclude_species_list",
+        ConfigValue(
+            default=None,
+            domain=list,
+            description="List of species to exclude from speciation in reaktoro",
+            doc="""The species that should not be included when speciating aqueous, liquid, or condensed phases, this 
+            will completely exclude the species from thermodynamic speciation. This is not the same as
+            rkt.exclude. """,
+        ),
+    )
+    CONFIG.declare(
+        "speciation_block_exclude_species",
+        ConfigValue(
+            default=None,
+            domain=list,
+            description="List of species to not include in output of speciation block and input into property block",
+            doc="""If provided then selected species will not be passed from speciation block and into property block. 
+            This is useful for removing species that do not go above 0 and cause issues during solve. No warning or checks
+            are performed by Reaktoro-PSE for scenario where the given species is above 0 at any time during solve. The user
+            should manually verify that the specie did not go above 0 or expected value. """,
+        ),
+    )
     CONFIG.declare(
         "exact_speciation",
         ConfigValue(
@@ -220,6 +247,19 @@ class ReaktoroBlockData(ProcessBlockData):
             """,
         ),
     )
+    CONFIG.declare(
+        "reaktoro_block_manager",
+        ConfigValue(
+            default=None,
+            domain=IsInstance(ReaktoroBlockManager),
+            description="Reaktoro block manager for parallelizing reaktoro solves",
+            doc="""
+            Option to provide a reaktoro block manager which would manage all blocks built on a model and
+            allow their parallel execution. When using ReaktorBlockManager, make sure to run
+            ReaktoroBlockManager.build_reaktoro_blocks() after constructing all of ReaktoroBlocks and 
+            before doing any initialization calls, or interacting with ReaktoroBlocks themself .""",
+        ),
+    )
     CONFIG.declare("jacobian_options", JacobianOptions().get_dict())
 
     CONFIG.declare(
@@ -327,7 +367,9 @@ def get_phases(phase_type):
                 building_prop_block_after_speciation()
                 and getattr(self.config, phase_type).phase_components is None
             ):
-                return getattr(self.speciation_block.rkt_state, phase_type)
+                return self.speciation_block.rkt_state.phase_manager.registered_phases[
+                    phase_type
+                ].phase_list
             else:
                 return getattr(self.config, phase_type).phase_components
 
@@ -368,8 +410,17 @@ def get_phases(phase_type):
         liquid_input_composition = self.config.liquid_phase.composition
         condensed_input_composition = self.config.condensed_phase.composition
         if building_prop_block_after_speciation():
+            # we need to ensure when we provide intial input compo into
+            # specitaiton block we don't have exteremely high ion concetration
+            # these value swill be overwritten during intilization anyway
+            for ion, obj in self.speciation_block.outputs.items():
+                if self.config.aqueous_phase.fixed_solvent_specie in ion:
+                    obj.set_value(obj.value * 10)
+                else:
+                    obj.set_value(obj.value / 1000)
             if aqueous_input_composition is not {}:
                 aqueous_input_composition = self.speciation_block.outputs
+
                 liquid_input_composition = {}
                 condensed_input_composition = {}
             elif liquid_input_composition is not {}:
@@ -386,6 +437,7 @@ def get_phases(phase_type):
                 )
         else:
             aqueous_input_composition = self.config.aqueous_phase.composition
+        block.rkt_state.register_species_to_exclude(self.config.exclude_species_list)
 
         block.rkt_state.register_aqueous_inputs(
             composition=aqueous_input_composition,
@@ -429,8 +481,12 @@ def get_phases(phase_type):
             block.rkt_state.register_aqueous_phase(get_phases("aqueous_phase"))
             block.rkt_state.register_liquid_phase(get_phases("liquid_phase"))
             block.rkt_state.register_solid_phases(get_phases("solid_phase"))
-            block.rkt_state.register_condensed_phase(get_phases("condensed_phase"))
-            block.rkt_state.register_gas_phase(get_phases("gas_phase"))
+            block.rkt_state.register_condensed_phase(
+                get_phases("condensed_phase"),
+            )
+            block.rkt_state.register_gas_phase(
+                get_phases("gas_phase"), self.config.gas_phase.speciate_phase_component
+            )
             block.rkt_state.register_mineral_phases(get_phases("mineral_phase"))
             block.rkt_state.register_ion_exchange_phase(
                 get_phases("ion_exchange_phase")
@@ -551,6 +607,7 @@ def build_rkt_inputs(
                 dissolve_species_in_rkt=self.config.dissolve_species_in_reaktoro,
                 exact_speciation=True,
             )
+        block.rkt_inputs.build_input_specs()
 
     def build_rkt_outputs(self, block, speciation_block=False):
         """this will build rkt outputs specified block.
@@ -570,7 +627,14 @@ def build_rkt_outputs(self, block, speciation_block=False):
             raise ValueError("Outputs must be provided!")
         if speciation_block:
             # when speciating we only want species amounts as output
-            block.rkt_outputs.register_output("speciesAmount", get_all_indexes=True)
+            if self.config.speciation_block_exclude_species is not None:
+                block.rkt_outputs.register_output(
+                    "speciesAmount",
+                    get_all_indexes=True,
+                    ignore_indexes=self.config.speciation_block_exclude_species,
+                )
+            else:
+                block.rkt_outputs.register_output("speciesAmount", get_all_indexes=True)
         else:
             # build user requested outputs
             for output_key, output_var in self.config.outputs.items():
@@ -584,7 +648,13 @@ def build_rkt_outputs(self, block, speciation_block=False):
                         output_prop = None
                     if isinstance(output_var, bool):
                         block.rkt_outputs.register_output(
-                            output_key, get_all_indexes=output_var
+                            output_key, get_all_indexes=True
+                        )
+                    elif isinstance(output_var, list):
+                        block.rkt_outputs.register_output(
+                            output_key,
+                            get_all_indexes=True,
+                            ignore_indexes=output_var,
                         )
                     else:
                         block.rkt_outputs.register_output(
@@ -663,7 +733,18 @@ def build_gray_box(self, block):
         block.rkt_block_builder.configure_jacobian_scaling(
             jacobian_scaling_type=scaling_type, user_scaling=scaling
         )
-        block.rkt_block_builder.build_reaktoro_block()
+        if self.config.reaktoro_block_manager is not None:
+            managed_block = self.config.reaktoro_block_manager.register_block(
+                state=block.rkt_state,
+                inputs=block.rkt_inputs,
+                outputs=block.rkt_outputs,
+                jacobian=block.rkt_jacobian,
+                solver=block.rkt_solver,
+                builder=block.rkt_block_builder,
+            )
+            block.managed_block = managed_block
+        else:
+            block.rkt_block_builder.build_reaktoro_block()
 
     # TODO: Update to provide output location (e.g. StringIO)
     def display_jacobian_outputs(self):
@@ -694,21 +775,24 @@ def display_reaktoro_state(self):
         """Displays reaktoro state"""
         if self.config.build_speciation_block:
             _log.info("-----Displaying information for speciation block ------")
-            _log.info(self.speciation_block.rkt_state.state)
+            self.speciation_block.rkt_block_builder.display_state()
         _log.info("-----Displaying information for property block ------")
-        _log.info(self.rkt_state.state)
+        self.rkt_block_builder.display_state()
 
-    def set_jacobian_scaling(self, user_scaling_dict, set_on_speciation_block=True):
-        """Sets jacobian scaling
+    def update_jacobian_scaling(
+        self, user_scaling_dict=None, set_on_speciation_block=True
+    ):
+        """This will recalculate jacobian scaling and update with user provided values
         Keywords:
         user_scaling_dict -- Dictionary that contains jacobian keys and scaling block
         set_on_speciation_block -- if scaling should be also set on speciation block if built.
         """
         if self.config.build_speciation_block and set_on_speciation_block:
+            self.speciation_block.rkt_block_builder.get_jacobian_scaling()
             self.speciation_block.rkt_block_builder.set_user_jacobian_scaling(
                 user_scaling_dict
             )
-
+        self.rkt_block_builder.get_jacobian_scaling()
         self.rkt_block_builder.set_user_jacobian_scaling(user_scaling_dict)
 
     # TODO: Update to use new initialization method https://idaes-pse.readthedocs.io/en/stable/reference_guides/initialization/developing_initializers.html?highlight=Initializer
diff --git a/src/reaktoro_pse/reaktoro_block_config/input_options.py b/src/reaktoro_pse/reaktoro_block_config/input_options.py
index 38f3ab3..1431b6f 100644
--- a/src/reaktoro_pse/reaktoro_block_config/input_options.py
+++ b/src/reaktoro_pse/reaktoro_block_config/input_options.py
@@ -7,7 +7,11 @@ def __init__(self):
         pass
 
     def get_dict(
-        self, include_pH=False, aqueous_phase=False, include_solvent_species=False
+        self,
+        include_pH=False,
+        aqueous_phase=False,
+        include_solvent_species=False,
+        include_speciate_phase_component=False,
     ):
         phase_input = ConfigDict()
         phase_input.declare(
@@ -44,6 +48,20 @@ def get_dict(
             """,
             ),
         )
+        if include_speciate_phase_component:
+            phase_input.declare(
+                "speciate_phase_component",
+                ConfigValue(
+                    default=False,
+                    description="To speciate specified phases",
+                    doc="""
+                        Use this to speciate give specie or elements to all possible
+                        phases or elements
+                        
+                        Works only on 
+                        Gas phase""",
+                ),
+            )
         phase_input.declare(
             "convert_to_rkt_species",
             ConfigValue(
diff --git a/src/reaktoro_pse/reaktoro_block_config/jacobian_options.py b/src/reaktoro_pse/reaktoro_block_config/jacobian_options.py
index 03074e4..f4a3890 100644
--- a/src/reaktoro_pse/reaktoro_block_config/jacobian_options.py
+++ b/src/reaktoro_pse/reaktoro_block_config/jacobian_options.py
@@ -30,7 +30,7 @@ def get_dict(self):
         CONFIG.declare(
             "numerical_order",
             ConfigValue(
-                default=8,
+                default=10,
                 domain=int,
                 description="Defines order of numerical jacobian (should be an even number)",
                 doc="""
@@ -43,7 +43,7 @@ def get_dict(self):
         CONFIG.declare(
             "numerical_step",
             ConfigValue(
-                default=1e-5,
+                default=1e-4,
                 domain=float,
                 description="Defines the step to use for numerical descritiazaiton",
                 doc="""This will define how small of a step to use for numerical derivative propagation which takes
@@ -89,7 +89,7 @@ def get_dict(self):
                 domain=IsInstance((str, HessTypes)),
                 description="Hessian type to use for reaktor gray box",
                 doc="""Hessian type to use, some might provide better stability
-                options (Jt.J, BFGS, BFGS-mod,BFGS-damp, BFGS-ipopt""",
+                options (Jt.J, BFGS, BFGS-mod, BFGS-damp, BFGS-ipopt""",
             ),
         )
         return CONFIG
diff --git a/src/reaktoro_pse/tests/test_reaktoro_block.py b/src/reaktoro_pse/tests/test_reaktoro_block.py
index d857460..a51984d 100644
--- a/src/reaktoro_pse/tests/test_reaktoro_block.py
+++ b/src/reaktoro_pse/tests/test_reaktoro_block.py
@@ -237,7 +237,12 @@ def test_blockBuild_with_speciation_block(build_rkt_state_with_species):
             pytest.approx(scaling_result["speciation_block"][key], 1e-3)
             == expected_scaling["speciation_block"][key]
         )
-    m.property_block.set_jacobian_scaling(new_scaling)
+    m.property_block.update_jacobian_scaling(new_scaling)
+    scaling_result = m.property_block.display_jacobian_scaling()
+
+    assert "speciation_block" in scaling_result
+    for key in scaling_result["speciation_block"]:
+        assert scaling_result["speciation_block"][key] == 1
     new_scaling = {}
     for key in scaling_result["property_block"]:
         new_scaling[key] = 1
@@ -245,12 +250,10 @@ def test_blockBuild_with_speciation_block(build_rkt_state_with_species):
             pytest.approx(scaling_result["property_block"][key], 1e-3)
             == expected_scaling["property_block"][key]
         )
-    m.property_block.set_jacobian_scaling(new_scaling)
+    m.property_block.update_jacobian_scaling(new_scaling)
     scaling_result = m.property_block.display_jacobian_scaling()
-    assert "speciation_block" in scaling_result
+
     assert "property_block" in scaling_result
-    for key in scaling_result["speciation_block"]:
-        assert scaling_result["speciation_block"][key] == 1
     for key in scaling_result["property_block"]:
         assert scaling_result["property_block"][key] == 1
 
diff --git a/src/reaktoro_pse/tutorials/basic_reaktoro_block_interaction.ipynb b/src/reaktoro_pse/tutorials/basic_reaktoro_block_interaction.ipynb
index 56ab8c5..de37e17 100644
--- a/src/reaktoro_pse/tutorials/basic_reaktoro_block_interaction.ipynb
+++ b/src/reaktoro_pse/tutorials/basic_reaktoro_block_interaction.ipynb
@@ -233,6 +233,7 @@
    "outputs": [],
    "source": [
     "\"\"\"Ions\"\"\"\n",
+    "\n",
     "ions = list(sea_water_composition.keys())\n",
     "\n",
     "\"\"\"Ions concentration variable\"\"\"\n",
@@ -259,6 +260,7 @@
    "outputs": [],
    "source": [
     "\"\"\"Charge neutrality\"\"\"\n",
+    "\n",
     "m.fs.sea_water.charge = Var(initialize=0)\n",
     "set_scaling_factor(m.fs.sea_water.charge, 1e8)\n",
     "\n",
@@ -362,7 +364,7 @@
     "        \"fixed_solvent_specie\": \"H2O\",  # We need to define our aqueous solvent as we have to speciate the block\n",
     "    },\n",
     "    outputs=m.fs.sea_water.outputs,  # outputs we desired\n",
-    "    database='PhreeqcDatabase', #can also be reaktoro.PhreeqcDatabase('pitzer.dat')\n",
+    "    database=\"PhreeqcDatabase\",  # can also be reaktoro.PhreeqcDatabase('pitzer.dat')\n",
     "    database_file=\"pitzer.dat\",  # needs to be a string that names the database file or points to its location\n",
     "    dissolve_species_in_reaktoro=True,  # This will sum up all species into elements in Reaktoro directly, if set to false, it will build Pyomo constraints instead\n",
     "    assert_charge_neutrality=False,  # This is True by Default, but here we actually want to adjust the input speciation till the charge is zero\n",
@@ -1065,7 +1067,7 @@
     "print(\"Sea water pH: \", m.fs.sea_water.pH.value)\n",
     "print(\"Reaktoro block outputs\")\n",
     "for key, obj in m.fs.modified_sea_water.outputs.items():\n",
-    "    print(key, obj.value) "
+    "    print(key, obj.value)"
    ]
   },
   {
@@ -1550,12 +1552,13 @@
    ],
    "source": [
     "\"\"\"Increase feed concentration\"\"\"\n",
+    "\n",
     "for ion, pyo_obj in m.fs.sea_water.species_concentrations.items():\n",
-    "    if ion!='Cl':\n",
-    "        pyo_obj.fix(pyo_obj.value*10)\n",
+    "    if ion != \"Cl\":\n",
+    "        pyo_obj.fix(pyo_obj.value * 10)\n",
     "\n",
     "m.fs.modified_sea_water.pH.unfix()\n",
-    "m.fs.modified_sea_water.outputs['scalingTendency','Calcite'].fix(1)\n",
+    "m.fs.modified_sea_water.outputs[\"scalingTendency\", \"Calcite\"].fix(1)\n",
     "m.fs.sea_water.species_concentrations.display()\n",
     "\n",
     "cy_solver = get_solver(solver=\"cyipopt-watertap\")\n",
diff --git a/src/reaktoro_pse/tutorials/integration_with_ro.ipynb b/src/reaktoro_pse/tutorials/integration_with_ro.ipynb
index 7d0a0b7..5efbb78 100644
--- a/src/reaktoro_pse/tutorials/integration_with_ro.ipynb
+++ b/src/reaktoro_pse/tutorials/integration_with_ro.ipynb
@@ -689,7 +689,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "m.fs.feed.reaktoro_properties.set_jacobian_scaling({(\"density\", None): 1000})"
+    "m.fs.feed.reaktoro_properties.update_jacobian_scaling({(\"density\", None): 1000})"
    ]
   },
   {