📚 Documentation

⚛️ Reaction Syatem

ReactionSystem

Bases: ThermoLinkDB, ReferenceManager

Class to represent a system of chemical reactions.

Source code in PyReactLab/docs/reactionsystem.py

class ReactionSystem(ThermoLinkDB, ReferenceManager):
    """Class to represent a system of chemical reactions."""

    # NOTE: class variables
    __system_name = None
    __reactions = None
    # primary analysis result
    __reaction_analysis = None
    __reaction_list: Dict[str, Reaction] = {}

    # reference plugin
    _references = {}

    # overall reaction phase
    overall_reaction_phase = None

    def __init__(
        self,
        system_name: str,
        reactions: List[Dict[str, Any]],
        model_source: Dict[str, Any],
        **kwargs
    ):
        self.__system_name = system_name
        self.__reactions = reactions

        # NOTE: model source
        self.__model_source = model_source

        # NOTE: kwargs
        # phase rule
        self.phase_rule = kwargs.get("phase_rule", None)

        # NOTE: init class
        ReferenceManager.__init__(self)
        ThermoLinkDB.__init__(self, model_source)

        # SECTION: init class
        self.ReactionAnalyzer_ = ReactionAnalyzer()

        # SECTION: load reference
        # reference plugin (default app params)
        self._references = self.load_reference()

        # SECTION: energy analysis result list
        self.__reaction_analyzer()

    @property
    def system_name(self) -> str:
        """Get the name of the reaction system."""
        # check
        if self.__system_name is None:
            return "No system name found."
        return self.__system_name

    @property
    def reactions(self) -> List[Dict[str, str]]:
        """Get the reactions of the reaction system."""
        # check
        if self.__reactions is None:
            raise ValueError("No reactions found.")
        return self.__reactions

    @property
    def reaction_list(self) -> Dict[str, Reaction]:
        """Get the list of reactions (object) in the reaction system."""
        # check
        if not self.__reaction_list:
            raise ValueError("No reactions found in the reaction system.")
        return self.__reaction_list

    def select_reaction(self, reaction_name: str) -> Reaction:
        """
        Select a reaction from the reaction system.

        Parameters
        ----------
        reaction_name : str
            Name of the reaction.

        Returns
        -------
        Reaction or None
            Reaction object if found, otherwise None.
        """
        try:
            # NOTE: check if reaction name is valid
            if reaction_name not in self.__reaction_list:
                raise ValueError(f"Invalid reaction name: {reaction_name}")

            # reaction object
            reaction = self.__reaction_list.get(reaction_name, None)
            # check if reaction is valid
            if reaction is None or reaction == 'None':
                raise ValueError(
                    f"Invalid reaction object for {reaction_name}")
            return reaction
        except Exception as e:
            raise Exception(
                f"Error in ReactionSystem.select_reaction(): {str(e)}") from e

    def __reaction_analyzer(self) -> None:
        """
        Execute the primary analysis for the reaction system.
        """
        try:
            # NOTE: initialize
            ChemReactUtils_ = ChemReactUtils()
            # ReactionAnalyzer_ = ReactionAnalyzer()

            # SECTION: reaction system analysis
            # analyze reaction
            reaction_res = {}

            # looping through each reaction
            for item in self.reactions:
                # NOTE: create reaction
                r_ = Reaction(
                    self.datasource,
                    self.equationsource,
                    item,
                    phase_rule=self.phase_rule,)

                # NOTE: analyze reaction
                _res = r_.reaction_analysis_result

                # name
                name: str = item['name']
                # update
                reaction_res[name] = _res
                # set
                self.__reaction_list[name] = r_

            # NOTE: reaction numbers
            # self.reaction_numbers = len(reaction_res)

            # SECTION: analyze overall reaction
            # ! to set consumed, produced, and intermediate species
            # NOTE: version 1
            # res_0 = ChemReactUtils_.analyze_overall_reactions(
            #     self.reactions)
            # NOTE: version 2
            res_0 = ChemReactUtils_.analyze_overall_reactions_v2(
                reaction_res)

            # SECTION: set component
            # NOTE: version 1
            # res_1 = ChemReactUtils_.define_component_id(
            #     reaction_res)
            # NOTE: version 2
            res_1 = ChemReactUtils_.define_component_id_v2(
                reaction_res)

            # extract
            # ? component_list: list of components
            # ? component_dict: dict of component: id
            # ? comp_list: list of dict of component stoichiometry
            # ? comp_coeff: list of list of component stoichiometry
            # ? component_state_list: list of component name and state
            (
                component_list,
                component_dict,
                comp_list,
                comp_coeff,
                component_state_list
            ) = res_1

            # set stoichiometry transpose
            comp_coeff_t = np.array(comp_coeff).T

            # SECTION: overall reaction phase
            # NOTE: set reaction states
            overall_reaction_phases = []
            for item in reaction_res:
                # name
                name = item

                # NOTE: reaction state
                _res = self.__reaction_list[name].reaction_phase

                # save
                overall_reaction_phases.append(_res)

            # set
            overall_reaction_phases = list(set(overall_reaction_phases))

            # check if all reactions have the same phase
            if len(overall_reaction_phases) == 0:
                raise ValueError(
                    "No overall reaction phases found in the reactions.")

            # set
            if len(overall_reaction_phases) == 1:
                # set
                self.overall_reaction_phase = overall_reaction_phases[0]
            else:
                # set
                self.overall_reaction_phase = '-'.join(overall_reaction_phases)

            # SECTION: phase components
            phase_contents = ChemReactUtils_.reaction_phase_analysis(
                reaction_res=reaction_res,
            )

            # SECTION: energy analysis
            # energy analysis result list
            energy_analysis = {}

            # loop through each reaction
            for item in reaction_res:
                # name
                name = item

                # NOTE: energy analysis
                reaction_ = self.__reaction_list[name]
                # ! energy analysis to get the following results:
                # ? gibbs free energy of formation (std)
                # ? enthalpy of formation (std)
                # ? gibbs free energy of reaction (std)
                # ? enthalpy of reaction (std)
                # ? equilibrium constant (std)
                _res = reaction_.energy_analysis

                # save
                energy_analysis[item] = _res

            # NOTE: set primary analysis result
            self.reaction_numbers = len(reaction_res)
            self.reaction_analysis = reaction_res
            self.overall_reaction_analysis = res_0
            self.component_list = component_list
            self.component_dict = component_dict
            self.coeff_list_dict = comp_list
            self.coeff_list_list = comp_coeff
            self.coeff_T_list_list = comp_coeff_t  # transpose
            self.component_state_list = component_state_list  # component state list
            self.energy_analysis = energy_analysis  # energy analysis result
            self.phase_contents = phase_contents  # phase contents

        except Exception as e:
            raise Exception(
                f"Error in ReactionSystem.go(): {str(e)}") from e

    def component_formation_energies(
        self,
        component_name: str,
        temperature: List[float | str],
        res_format: Literal[
            'symbolic', 'names'
        ] = 'names',
        message: Optional[str] = None
    ):
        '''
        Calculate enthalpy and gibbs free energy of formation for a component at a given temperature.

        Parameters
        ----------
        component_name : str
            Name of the component such as 'H2O-l', 'CO2-g', etc.
        temperature : list[float, str]
            Temperature in any unit, e.g. [300.0, "K"], It is automatically converted to Kelvin.
        res_format : str, optional
            Format of the result, by default 'names'.
        message : str, optional
            Optional message to display during the calculation, by default None.


        '''
        try:
            # NOTE: check if T
            # check if T is a list
            if not isinstance(temperature, list):
                raise ValueError("Temperature must be a list.")

            # check if T is a number
            if not isinstance(temperature[0], (int, float)):
                raise ValueError("Temperature must be a number.")

            # check if T is a string
            if not isinstance(temperature[1], str):
                raise ValueError("Temperature unit must be a string.")

            # NOTE: convert temperature to Kelvin
            # get
            T_value = temperature[0]
            T_unit = temperature[1]
            # set unit
            unit_set = f"{T_unit} => K"
            T = pycuc.to(T_value, unit_set)

            # SECTION: check if component name is valid
            # init
            component_name_ = []
            component_name = component_name.strip()

            # loop through component state list
            for item in self.component_state_list:
                # check if component name is valid
                if item[2] == component_name.strip():
                    # set
                    component_name_.append(item[0])
                    component_name_.append(item[2])  # state
                    break

            # ! check if component name is valid
            if len(component_name_) == 0:
                raise ValueError(
                    f"Invalid component name: {component_name}. Please check the component name and try again.")

            # SECTION: calculate component formation energies
            res_ = self.ReactionAnalyzer_.component_energy_at_temperature(
                datasource=self.datasource,
                equationsource=self.equationsource,
                component_names=component_name_,
                temperature=T,
                res_format=res_format,
            )

            # NOTE: add message if provided
            if message is None:
                message = f"Formation energies of {component_name_} at {T_value} {T_unit}"

            # set message
            res_['message'] = {'value': message}

            # return res
            return res_
        except Exception as e:
            raise Exception(
                f"Error in ReactionSystem.calc_component_formation_energies(): {str(e)}") from e

    def reaction_equilibrium_constant(
        self,
        reaction_name: str,
        temperature: list[float | str],
        method: Literal[
            "van't Hoff", "shortcut van't Hoff"
        ] = "van't Hoff",
        **kwargs
    ):
        """
        Calculate the equilibrium constant at a given temperature using the van't Hoff equation.

        Parameters
        ----------
        reaction_name : str
            Name of the reaction.
        temperature : list
            Temperature in the form of [value, unit], the unit is automatically converted to K.
        method : str, optional
            Method to calculate the equilibrium constant, by default "van't Hoff".
            Options are "van't Hoff" or "shortcut van't Hoff".
        **kwargs : dict
            Additional arguments for the calculation.
                - message: Optional message to display during the calculation.

        Returns
        -------
        res : dict
            Equilibrium constant at the given temperature.
        """
        try:
            # NOTE: check if temperature is valid
            if not isinstance(temperature, list):
                raise ValueError("Temperature must be a number.")

            # check if temperature is valid
            if len(temperature) != 2:
                raise ValueError(
                    "Temperature must be a list of length 2 following [value, unit].")

            if not isinstance(temperature[0], (int, float)):
                raise ValueError("Temperature must be a number.")

            if not isinstance(temperature[1], str):
                raise ValueError("Temperature unit must be a string.")

            # NOTE: check if reaction name is valid
            if reaction_name not in self.__reaction_list:
                raise ValueError(f"Invalid reaction name: {reaction_name}")

            # SECTION: pre-calculation
            # get reaction object
            reaction = self.__reaction_list[reaction_name]

            # check if reaction is valid
            if not isinstance(reaction, Reaction):
                raise ValueError(
                    f"Invalid reaction object for {reaction_name}")

            # NOTE: calculate equilibrium constant at the given temperature
            res = reaction.cal_equilibrium_constant(
                temperature,
                method=method
            )

            # NOTE: message
            message = kwargs.get("message", None)
            if message is not None:
                res['message'] = message

            # res
            return res
        except Exception as e:
            raise Exception(
                f"Error in ReactionSystem.equilibrium_constant_at_temperature(): {str(e)}") from e

    def equilibrium(
        self,
        inputs: Dict[str, Any],
        conversion: Optional[List[str]] = None,
        gas_mixture: Literal[
            "ideal", "non-ideal"
        ] = "ideal",
        solution: Literal[
            "ideal", "non-ideal"
        ] = "ideal",
        method: Literal[
            'minimize', 'least_squares'
        ] = 'minimize',
        **kwargs
    ):
        """
        Calculate the equilibrium state of the reaction system.

        Parameters
        ----------
        inputs : dict
            Inputs for the equilibrium calculation.
        conversion : list, optional
            List of components to calculate conversion for, by default None.
        gas_mixture : str, optional
            Type of gas mixture, by default "ideal".
        solution : str, optional
            Type of liquid mixture, by default "ideal".
        method : str, optional
            Method for the calculation, by default "minimize".
            Options are "minimize" or "least_squares".
        **kwargs : dict
            Additional arguments for the calculation.
                - eos_model: Equation of state model to use for the calculation. Options are "SRK" or "PR".
                - activity_model: Activity model to use for the calculation. Options are "NRTL" or "UNIFAC".
                - message: Optional message to display during the calculation.


        Returns
        -------
        dict
            Equilibrium state of the reaction system.

        Notes
        -----
        Inputs must contain `mole_fraction` or `mole`, `temperature`, and `pressure` as follows:

        - mole_fraction: dict, initial mole fraction of the components in the system.
        - mole: dict, initial mole of the components in the system.
        - temperature: list, temperature in the form of [value, unit].
        - pressure: list, pressure in the form of [value, unit].

        For `non-ideal solution`, the rules of the activity model should be followed.

        Activity inputs for NRTL model should contain:
        - tau: dict, tau values for the components.
        - alpha: dict, alpha values for the components.

        Activity inputs for UNIQUAC model should contain:
        - q: dict, q values for the components.
        - r: dict, r values for the components.
        - tau: dict, tau values for the components.

        For both models, dU or a, b, c, and d are required to calculate tau in case tau is not provided. All values should be introduced in the `inputs` as:

        ```python
        inputs = {
            'mole': mole,
            'temperature': [100, "C"],
            'pressure': [1.0, "bar"],
            'tau': tau,
            'alpha': alpha, ...
            }
        ```
        """
        try:
            # ! Start timing
            start_time = time.time()

            # SECTION: check args
            # NOTE: check if inputs are valid
            if not isinstance(inputs, dict):
                raise ValueError("Inputs must be a dictionary.")
            # check if inputs are valid
            if not any(key in inputs for key in ["mole_fraction", "mole"]):
                raise ValueError(
                    "Inputs must contain either mole_fraction or mole.")
            if "temperature" not in inputs:
                raise ValueError("Inputs must contain temperature.")
            if "pressure" not in inputs:
                raise ValueError("Inputs must contain pressure.")

            # NOTE: check if temperature is valid
            # set
            temperature = inputs["temperature"]
            # check if temperature is valid
            if not isinstance(temperature, list):
                raise ValueError("Temperature must be a number.")

            # check if temperature is valid
            if len(temperature) != 2:
                raise ValueError(
                    "Temperature must be a list of length 2 following [value, unit].")

            if not isinstance(temperature[0], (int, float)):
                raise ValueError("Temperature must be a number.")

            if not isinstance(temperature[1], str):
                raise ValueError("Temperature unit must be a string.")

            # ! convert to **K**
            # set unit
            unit_set = f"{temperature[1]} => K"
            temperature_K = pycuc.to(temperature[0], unit_set)

            # NOTE: check if pressure is valid
            # set
            pressure = inputs["pressure"]
            # check if pressure is valid
            if not isinstance(pressure, list):
                raise ValueError("Pressure must be a number.")

            # check if pressure is valid
            if len(pressure) != 2:
                raise ValueError(
                    "Pressure must be a list of length 2 following [value, unit].")

            if not isinstance(pressure[0], (int, float)):
                raise ValueError("Pressure must be a number.")

            if not isinstance(pressure[1], str):
                raise ValueError("Pressure unit must be a string.")

            # ! convert to **bar**
            # set unit
            unit_set = f"{pressure[1]} => bar"
            pressure_bar = pycuc.to(pressure[0], unit_set)

            # NOTE: check if initial mole fraction is valid
            # set
            initial_mole_fraction = inputs.get("mole_fraction", None)
            # check
            initial_mole = inputs.get("mole", None)

            if initial_mole_fraction is not None:
                # check if initial mole fraction is valid
                if not isinstance(initial_mole_fraction, dict):
                    raise ValueError(
                        "Initial mole fraction must be a dictionary.")

                # check if initial mole fraction is valid
                for key in initial_mole_fraction:
                    if not isinstance(key, str):
                        raise ValueError(
                            "Initial mole fraction key must be a string.")
                    if not isinstance(initial_mole_fraction[key], (int, float)):
                        raise ValueError(
                            "Initial mole fraction value must be a number.")

                # ? normalize
                # check if sum of mole fraction is 1
                initial_mole_fraction = ReactionAnalyzer.norm_mole_fraction(
                    initial_mole_fraction)

                # NOTE: convert to mole fraction
                initial_mole = ReactionAnalyzer.cal_mole(
                    initial_mole_fraction)

            elif initial_mole is not None:
                # check if initial mole is valid
                if not isinstance(initial_mole, dict):
                    raise ValueError("Initial mole must be a dictionary.")

                # check if initial mole is valid
                for key in initial_mole:
                    if not isinstance(key, str):
                        raise ValueError(
                            "Initial mole key must be a string.")
                    if not isinstance(initial_mole[key], (int, float)):
                        raise ValueError(
                            "Initial mole value must be a number.")

                # NOTE: convert to mole fraction
                initial_mole_fraction, _ = ReactionAnalyzer.cal_mole_fraction(
                    initial_mole)

            # NOTE: build mole and mole fraction matrix regarding the component dict
            # check
            if (initial_mole is not None and
                    initial_mole_fraction is not None):
                # set values
                (
                    initial_mole_std, initial_mole_fraction_std, _, _
                ) = ReactionAnalyzer.set_stream(
                    component_dict=self.component_dict,
                    mole=initial_mole,
                    mole_fraction=initial_mole_fraction,
                )
            else:
                raise ValueError(
                    "Initial mole and mole fraction must be provided.")

            # SECTION: kwargs
            # eos model (name)
            eos_model = kwargs.get("eos_model", "SRK")
            # check if eos model is valid
            if eos_model not in ["SRK", "PR", 'RK', 'vdW']:
                raise ValueError(
                    "Invalid eos model. Options are 'SRK','RK','PR', or 'vdW'.")

            # activity model (name)
            activity_model = kwargs.get("activity_model", "NRTL")
            # check if activity model is valid
            if activity_model not in ["NRTL", "UNIQUAC"]:
                raise ValueError(
                    "Invalid activity model. Options are 'NRTL' or 'UNIQUAC'.")

            # SECTION: init
            ReactionOptimizer_ = ReactionOptimizer(
                self.datasource,
                self.equationsource,
                self.component_dict,
                self.coeff_list_dict,
                self.coeff_T_list_list,
                self.reaction_analysis,
                self.overall_reaction_analysis,
            )

            # SECTION: setting up the reaction optimizer
            # NOTE: set up eos model
            # ? eos model
            ReactionOptimizer_.eos_model = eos_model
            # init eos class
            ReactionOptimizer_.eos = ptm.eos()

            # gas mixture
            ReactionOptimizer_.gas_mixture = gas_mixture

            # NOTE: set up activity model
            # ? activity inputs
            # NRTL: tau, alpha,
            # dg or a, b, c, and d are required to calculated tau
            # UNIQUAC: q, r, tau
            # dU or a, b, c, and d are required to calculated tau
            activity_inputs = {**inputs}
            # remove temperature and pressure
            activity_inputs.pop("temperature", None)
            activity_inputs.pop("pressure", None)
            # remove mole and mole fraction
            activity_inputs.pop("mole", None)
            activity_inputs.pop("mole_fraction", None)

            # set
            if activity_inputs is not None:
                #  check if activity inputs is valid
                if not isinstance(activity_inputs, dict):
                    raise ValueError(
                        "Activity inputs must be a dictionary.")

                # set
                ReactionOptimizer_.activity_inputs = activity_inputs

            # ? activity model
            ReactionOptimizer_.activity_model = activity_model
            # init activity class
            ReactionOptimizer_.activity = ptm.activities(
                components=self.component_list,
                model_name=activity_model,
                model_source=self.model_source,
            )

            # solution
            ReactionOptimizer_.solution = solution

            # SECTION: equilibrium constant calculation
            # res
            equilibrium_constant = {}

            # loop through each reaction
            for key, r in self.__reaction_list.items():
                # NOTE: check if reaction is valid
                if not isinstance(r, Reaction):
                    raise ValueError(
                        f"Invalid reaction object for {key}")

                # ! calculate equilibrium constant at the given temperature
                res_ = r.cal_equilibrium_constant(temperature)

                # save
                equilibrium_constant[key] = res_

            # NOTE: run equilibrium calculation
            opt_res = ReactionOptimizer_.opt_run(
                initial_mole=initial_mole_std,
                initial_mole_fraction=initial_mole_fraction_std,
                temperature=temperature_K,
                pressure=pressure_bar,
                equilibrium_constants=equilibrium_constant,
                reaction_numbers=self.reaction_numbers,
                method=method,
            )

            # NOTE: check optimization result
            res = {}

            # check if optimization is successful
            if opt_res.success:
                # extent of reaction
                EoR = opt_res.x

                # calculate equilibrium state
                Eq_Xfs, Eq_Xfs_vector, Eq_Nf, Eq_Nfs_vector, Eq_Nfs = \
                    ReactionOptimizer_.equilibrium_results(
                        initial_mole=initial_mole_std,
                        EoR=EoR,
                    )

                # NOTE: calculate conversion
                # res
                conversion_res = None
                # conversion
                if conversion is not None:
                    # check list
                    if not isinstance(conversion, list):
                        raise ValueError(
                            "Conversion must be a list of strings.")

                    # check component in the components
                    for key in conversion:
                        if key not in self.component_dict:
                            raise ValueError(
                                f"Invalid component: {key} in the components.")

                    # calculate conversion
                    conversion_res = ReactionAnalyzer.cal_conversion(
                        initial_mole=initial_mole_std,
                        final_mole=Eq_Nfs,
                        components=conversion,
                    )

                # SECTION: set results
                # initial feed
                res['feed'] = {
                    "mole": {
                        'value': initial_mole_std,
                        'unit': "mol"
                    },
                    "mole_total": {
                        'value': sum(initial_mole_std.values()),
                        'unit': "mol"
                    },
                    "mole_fraction": {
                        'value': initial_mole_fraction_std,
                        'unit': "dimensionless"
                    },
                    "mole_fraction_sum": {
                        'value': sum(initial_mole_fraction_std.values()),
                        'unit': "dimensionless"
                    },
                }

                # equilibrium state
                res['equilibrium'] = {
                    "mole": {
                        'value': Eq_Nfs,
                        'unit': "mol"
                    },
                    "mole_total": {
                        'value': sum(Eq_Nfs.values()),
                        'unit': "mol"
                    },
                    "mole_fraction": {
                        'value': Eq_Xfs,
                        'unit': "dimensionless"
                    },
                    "mole_fraction_sum": {
                        'value': sum(Eq_Xfs.values()),
                        'unit': "dimensionless"
                    },
                }

                # extent of reaction
                res['extent_of_reaction'] = {
                    "value": EoR,
                    "unit": "dimensionless"
                }

                # equilibrium condition
                # temperature
                res['temperature'] = {
                    "value": temperature_K,
                    "unit": "K"
                }

                # pressure
                res['pressure'] = {
                    "value": pressure_bar,
                    "unit": "bar"
                }

                # optimization fun value
                res['optimization_fun'] = {
                    "value": opt_res.fun,
                    "unit": "dimensionless"
                }

                # conversion
                if conversion_res is not None:
                    res['conversion'] = conversion_res

                # message
                message = kwargs.get("message", None)
                if message is not None:
                    res['message'] = message

                # gas mixture
                res['gas_mixture'] = gas_mixture
                # solution
                res['solution'] = solution

            # NOTE: set time
            # ! Stop timing
            end_time = time.time()
            computation_time = end_time - start_time
            # add to res
            res['computation_time'] = {
                "value": computation_time,
                "unit": "s"
            }

            # res
            return res

        except Exception as e:
            raise Exception(
                f"Failing in the equilibrium calculations {str(e)}") from e

    def equilibrium_deviation(
        self,
        inputs: Dict[str, Any],
        gas_mixture: Literal[
            "ideal", "non-ideal"
        ] = "ideal",
        solution: Literal[
            "ideal", "non-ideal"
        ] = "ideal",
        **kwargs
    ):
        """
        Conduct equilibrium deviation calculations for the reaction system.

        This method calculates `the actual Gibbs free energy of reaction` (GiEn_rxn_T) at a given temperature and pressure, and component composition.

        There are three cases as:
        - Case 1: `GiEn_rxn_T < 0`, the reaction is not at equilibrium and the reaction will proceed in the forward direction.
        - Case 2: `GiEn_rxn_T > 0`, the reaction is not at equilibrium and the reaction will proceed in the reverse direction.
        - Case 3: `GiEn_rxn_T = 0`, the reaction is at equilibrium and no further reaction will occur.

        Parameters
        ----------
        inputs : dict
            Inputs for the equilibrium deviation calculation.
        conversion : list, optional
            List of components to calculate conversion for, by default None.
        gas_mixture : str, optional
            Type of gas mixture, by default "ideal".
        solution : str, optional
            Type of liquid mixture, by default "ideal".
        **kwargs : dict
            Additional arguments for the calculation.
                - eos_model: Equation of state model to use for the calculation. Options are "SRK" or "PR".
                - activity_model: Activity model to use for the calculation. Options are "NRTL" or "UNIFAC".
                - message: Optional message to display during the calculation.
                - mole_basis: 1 mole
                - minimum_mole: Minimum mole for the initial mole, by default 1e-5.

        Returns
        -------
        dict
            Equilibrium deviation of the reaction system.

        Notes
        -----
        The inputs and kwargs are similar to those in the `equilibrium` method.
        """
        # NOTE: call equilibrium method with deviation=True
        try:
            # ! Start timing
            start_time = time.time()

            if not self.overall_reaction_phase:
                raise ValueError(
                    "Overall reaction phase is not set. Please run the primary analysis first.")

            # mole basis
            mole_basis = kwargs.get("mole_basis", 1.0)
            # minimum mole
            minimum_mole = kwargs.get("minimum_mole", 1e-5)

            # SECTION: check args
            # NOTE: check if inputs are valid
            if not isinstance(inputs, dict):
                raise ValueError("Inputs must be a dictionary.")
            # check if inputs are valid
            if not any(key in inputs for key in ["mole_fraction", "mole"]):
                raise ValueError(
                    "Inputs must contain either mole_fraction or mole.")
            if "temperature" not in inputs:
                raise ValueError("Inputs must contain temperature.")
            if "pressure" not in inputs:
                raise ValueError("Inputs must contain pressure.")

            # SECTION: set temperature and pressure
            # NOTE: check if temperature is valid
            # set
            temperature = inputs["temperature"]
            # check if temperature is valid
            if not isinstance(temperature, list):
                raise ValueError("Temperature must be a number.")

            # check if temperature is valid
            if len(temperature) != 2:
                raise ValueError(
                    "Temperature must be a list of length 2 following [value, unit].")

            if not isinstance(temperature[0], (int, float)):
                raise ValueError("Temperature must be a number.")

            if not isinstance(temperature[1], str):
                raise ValueError("Temperature unit must be a string.")

            # ! convert to **K**
            # set unit
            unit_set = f"{temperature[1]} => K"
            temperature_K = pycuc.to(temperature[0], unit_set)
            # set temperature
            temperature = Temperature(
                value=temperature_K,
                unit="K"
            )

            # NOTE: check if pressure is valid
            # set
            pressure = inputs["pressure"]
            # check if pressure is valid
            if not isinstance(pressure, list):
                raise ValueError("Pressure must be a number.")

            # check if pressure is valid
            if len(pressure) != 2:
                raise ValueError(
                    "Pressure must be a list of length 2 following [value, unit].")

            if not isinstance(pressure[0], (int, float)):
                raise ValueError("Pressure must be a number.")

            if not isinstance(pressure[1], str):
                raise ValueError("Pressure unit must be a string.")

            # ! convert to **bar**
            # set unit
            unit_set = f"{pressure[1]} => bar"
            pressure_bar = pycuc.to(pressure[0], unit_set)
            # set pressure
            pressure = Pressure(
                value=pressure_bar,
                unit="bar"
            )

            # NOTE: set operating conditions
            operating_conditions = OperatingConditions(
                temperature=temperature,
                pressure=pressure,
            )

            # SECTION:
            # NOTE: check if initial mole fraction is valid
            # set
            initial_mole_fraction = inputs.get("mole_fraction", None)
            # check
            initial_mole = inputs.get("mole", None)

            if initial_mole_fraction is not None:
                # check if initial mole fraction is valid
                if not isinstance(initial_mole_fraction, dict):
                    raise ValueError(
                        "Initial mole fraction must be a dictionary.")

                # check if initial mole fraction is valid
                for key in initial_mole_fraction:
                    if not isinstance(key, str):
                        raise ValueError(
                            "Initial mole fraction key must be a string.")
                    if not isinstance(initial_mole_fraction[key], (int, float)):
                        raise ValueError(
                            "Initial mole fraction value must be a number.")

                # check if sum of mole fraction is 1
                initial_mole_fraction = ReactionAnalyzer.norm_mole_fraction(
                    initial_mole_fraction)

                # NOTE: convert to mole fraction
                initial_mole = ReactionAnalyzer.cal_mole(
                    initial_mole_fraction,
                    mole_basis
                )

            elif initial_mole is not None:
                # check if initial mole is valid
                if not isinstance(initial_mole, dict):
                    raise ValueError("Initial mole must be a dictionary.")

                # check if initial mole is valid
                for key in initial_mole:
                    if not isinstance(key, str):
                        raise ValueError(
                            "Initial mole key must be a string.")
                    if not isinstance(initial_mole[key], (int, float)):
                        raise ValueError(
                            "Initial mole value must be a number.")

                # set initial mole
                initial_mole = ReactionAnalyzer.set_initial_mole(
                    initial_mole,
                    minimum_mole=minimum_mole
                )

                # NOTE: convert to mole fraction
                initial_mole_fraction, _ = ReactionAnalyzer.cal_mole_fraction(
                    initial_mole)

            # NOTE: build mole and mole fraction matrix regarding the component dict
            # check
            if initial_mole is not None and initial_mole_fraction is not None:
                # set values
                (
                    initial_mole_std, initial_mole_fraction_std, _, _
                ) = ReactionAnalyzer.set_stream(
                    component_dict=self.component_dict,
                    mole=initial_mole,
                    mole_fraction=initial_mole_fraction,
                )
            else:
                raise ValueError(
                    "Initial mole and mole fraction must be provided.")

            # SECTION: mole and mole fraction for gas mixture and solution
            self.phase_stream = ReactionAnalyzer.set_phase_stream(
                initial_mole=initial_mole_std,
                initial_mole_fraction=initial_mole_fraction_std,
                phase_contents=self.phase_contents,
            )

            # SECTION: kwargs
            # eos model (name)
            eos_model = kwargs.get("eos_model", "SRK")
            # check if eos model is valid
            if eos_model not in ["SRK", "PR", 'RK', 'vdW']:
                raise ValueError(
                    "Invalid eos model. Options are 'SRK','RK','PR', or 'vdW'.")

            # activity model (name)
            activity_model = kwargs.get("activity_model", "NRTL")
            # check if activity model is valid
            if activity_model not in ["NRTL", "UNIQUAC"]:
                raise ValueError(
                    "Invalid activity model. Options are 'NRTL' or 'UNIQUAC'.")

            # SECTION: init chemical potential
            # init
            ChemicalPotential_ = ChemicalPotential(
                self.datasource,
                self.equationsource,
                self.reaction_list,
                self.component_dict,
                self.component_state_list,
                self.coeff_list_dict,
                self.reaction_analysis,
                self.phase_stream,
                self.phase_contents,
                self.overall_reaction_analysis,
                self.overall_reaction_phase,
                operating_conditions=operating_conditions,
            )

            # SECTION: setting up the chemical potential
            # NOTE: set up eos model
            # ? eos model
            ChemicalPotential_.eos_model = eos_model
            # init eos class
            ChemicalPotential_.eos = ptm.eos()

            # gas mixture
            ChemicalPotential_.gas_mixture = gas_mixture

            # NOTE: set up activity model
            # ? activity inputs
            # NRTL: tau, alpha,
            # dg or a, b, c, and d are required to calculated tau
            # UNIQUAC: q, r, tau
            # dU or a, b, c, and d are required to calculated tau
            activity_inputs = {**inputs}
            # remove temperature and pressure
            activity_inputs.pop("temperature", None)
            activity_inputs.pop("pressure", None)
            # remove mole and mole fraction
            activity_inputs.pop("mole", None)
            activity_inputs.pop("mole_fraction", None)

            # set
            if activity_inputs is not None:
                #  check if activity inputs is valid
                if not isinstance(activity_inputs, dict):
                    raise ValueError(
                        "Activity inputs must be a dictionary.")

                # set
                ChemicalPotential_.activity_inputs = activity_inputs

            # ? activity model
            ChemicalPotential_.activity_model = activity_model
            # ? set component list (liquid phase)
            # components_l = self.phase_contents['l']

            # init activity class
            ChemicalPotential_.activity = ptm.activities(
                components=self.component_list,
                model_name=activity_model,
                model_source=self.model_source,
            )

            # solution
            ChemicalPotential_.solution = solution

            # SECTION: the actual Gibbs free energy of reaction
            # calc
            res = ChemicalPotential_.cal_actual_gibbs_energy_of_reaction(
                operating_conditions=operating_conditions,
                gas_mixture=gas_mixture,
                solution=solution,
            )

            # NOTE: set time
            # ! Stop timing
            end_time = time.time()
            computation_time = end_time - start_time
            # add to res
            res['computation_time'] = {
                "value": computation_time,
                "unit": "s"
            }

            # res
            return res
        except Exception as e:
            raise Exception(
                f"Failing in the equilibrium deviation calculations {str(e)}") from e

reaction_list: Dict[str, Reaction] property

Get the list of reactions (object) in the reaction system.

reactions: List[Dict[str, str]] property

Get the reactions of the reaction system.

system_name: str property

Get the name of the reaction system.

__reaction_analyzer()

Execute the primary analysis for the reaction system.

Source code in PyReactLab/docs/reactionsystem.py

def __reaction_analyzer(self) -> None:
    """
    Execute the primary analysis for the reaction system.
    """
    try:
        # NOTE: initialize
        ChemReactUtils_ = ChemReactUtils()
        # ReactionAnalyzer_ = ReactionAnalyzer()

        # SECTION: reaction system analysis
        # analyze reaction
        reaction_res = {}

        # looping through each reaction
        for item in self.reactions:
            # NOTE: create reaction
            r_ = Reaction(
                self.datasource,
                self.equationsource,
                item,
                phase_rule=self.phase_rule,)

            # NOTE: analyze reaction
            _res = r_.reaction_analysis_result

            # name
            name: str = item['name']
            # update
            reaction_res[name] = _res
            # set
            self.__reaction_list[name] = r_

        # NOTE: reaction numbers
        # self.reaction_numbers = len(reaction_res)

        # SECTION: analyze overall reaction
        # ! to set consumed, produced, and intermediate species
        # NOTE: version 1
        # res_0 = ChemReactUtils_.analyze_overall_reactions(
        #     self.reactions)
        # NOTE: version 2
        res_0 = ChemReactUtils_.analyze_overall_reactions_v2(
            reaction_res)

        # SECTION: set component
        # NOTE: version 1
        # res_1 = ChemReactUtils_.define_component_id(
        #     reaction_res)
        # NOTE: version 2
        res_1 = ChemReactUtils_.define_component_id_v2(
            reaction_res)

        # extract
        # ? component_list: list of components
        # ? component_dict: dict of component: id
        # ? comp_list: list of dict of component stoichiometry
        # ? comp_coeff: list of list of component stoichiometry
        # ? component_state_list: list of component name and state
        (
            component_list,
            component_dict,
            comp_list,
            comp_coeff,
            component_state_list
        ) = res_1

        # set stoichiometry transpose
        comp_coeff_t = np.array(comp_coeff).T

        # SECTION: overall reaction phase
        # NOTE: set reaction states
        overall_reaction_phases = []
        for item in reaction_res:
            # name
            name = item

            # NOTE: reaction state
            _res = self.__reaction_list[name].reaction_phase

            # save
            overall_reaction_phases.append(_res)

        # set
        overall_reaction_phases = list(set(overall_reaction_phases))

        # check if all reactions have the same phase
        if len(overall_reaction_phases) == 0:
            raise ValueError(
                "No overall reaction phases found in the reactions.")

        # set
        if len(overall_reaction_phases) == 1:
            # set
            self.overall_reaction_phase = overall_reaction_phases[0]
        else:
            # set
            self.overall_reaction_phase = '-'.join(overall_reaction_phases)

        # SECTION: phase components
        phase_contents = ChemReactUtils_.reaction_phase_analysis(
            reaction_res=reaction_res,
        )

        # SECTION: energy analysis
        # energy analysis result list
        energy_analysis = {}

        # loop through each reaction
        for item in reaction_res:
            # name
            name = item

            # NOTE: energy analysis
            reaction_ = self.__reaction_list[name]
            # ! energy analysis to get the following results:
            # ? gibbs free energy of formation (std)
            # ? enthalpy of formation (std)
            # ? gibbs free energy of reaction (std)
            # ? enthalpy of reaction (std)
            # ? equilibrium constant (std)
            _res = reaction_.energy_analysis

            # save
            energy_analysis[item] = _res

        # NOTE: set primary analysis result
        self.reaction_numbers = len(reaction_res)
        self.reaction_analysis = reaction_res
        self.overall_reaction_analysis = res_0
        self.component_list = component_list
        self.component_dict = component_dict
        self.coeff_list_dict = comp_list
        self.coeff_list_list = comp_coeff
        self.coeff_T_list_list = comp_coeff_t  # transpose
        self.component_state_list = component_state_list  # component state list
        self.energy_analysis = energy_analysis  # energy analysis result
        self.phase_contents = phase_contents  # phase contents

    except Exception as e:
        raise Exception(
            f"Error in ReactionSystem.go(): {str(e)}") from e

component_formation_energies(component_name, temperature, res_format='names', message=None)

Calculate enthalpy and gibbs free energy of formation for a component at a given temperature.

Parameters

component_name : str Name of the component such as 'H2O-l', 'CO2-g', etc. temperature : list[float, str] Temperature in any unit, e.g. [300.0, "K"], It is automatically converted to Kelvin. res_format : str, optional Format of the result, by default 'names'. message : str, optional Optional message to display during the calculation, by default None.

Source code in PyReactLab/docs/reactionsystem.py

def component_formation_energies(
    self,
    component_name: str,
    temperature: List[float | str],
    res_format: Literal[
        'symbolic', 'names'
    ] = 'names',
    message: Optional[str] = None
):
    '''
    Calculate enthalpy and gibbs free energy of formation for a component at a given temperature.

    Parameters
    ----------
    component_name : str
        Name of the component such as 'H2O-l', 'CO2-g', etc.
    temperature : list[float, str]
        Temperature in any unit, e.g. [300.0, "K"], It is automatically converted to Kelvin.
    res_format : str, optional
        Format of the result, by default 'names'.
    message : str, optional
        Optional message to display during the calculation, by default None.


    '''
    try:
        # NOTE: check if T
        # check if T is a list
        if not isinstance(temperature, list):
            raise ValueError("Temperature must be a list.")

        # check if T is a number
        if not isinstance(temperature[0], (int, float)):
            raise ValueError("Temperature must be a number.")

        # check if T is a string
        if not isinstance(temperature[1], str):
            raise ValueError("Temperature unit must be a string.")

        # NOTE: convert temperature to Kelvin
        # get
        T_value = temperature[0]
        T_unit = temperature[1]
        # set unit
        unit_set = f"{T_unit} => K"
        T = pycuc.to(T_value, unit_set)

        # SECTION: check if component name is valid
        # init
        component_name_ = []
        component_name = component_name.strip()

        # loop through component state list
        for item in self.component_state_list:
            # check if component name is valid
            if item[2] == component_name.strip():
                # set
                component_name_.append(item[0])
                component_name_.append(item[2])  # state
                break

        # ! check if component name is valid
        if len(component_name_) == 0:
            raise ValueError(
                f"Invalid component name: {component_name}. Please check the component name and try again.")

        # SECTION: calculate component formation energies
        res_ = self.ReactionAnalyzer_.component_energy_at_temperature(
            datasource=self.datasource,
            equationsource=self.equationsource,
            component_names=component_name_,
            temperature=T,
            res_format=res_format,
        )

        # NOTE: add message if provided
        if message is None:
            message = f"Formation energies of {component_name_} at {T_value} {T_unit}"

        # set message
        res_['message'] = {'value': message}

        # return res
        return res_
    except Exception as e:
        raise Exception(
            f"Error in ReactionSystem.calc_component_formation_energies(): {str(e)}") from e

equilibrium(inputs, conversion=None, gas_mixture='ideal', solution='ideal', method='minimize', **kwargs)

Calculate the equilibrium state of the reaction system.

Parameters

inputs : dict Inputs for the equilibrium calculation. conversion : list, optional List of components to calculate conversion for, by default None. gas_mixture : str, optional Type of gas mixture, by default "ideal". solution : str, optional Type of liquid mixture, by default "ideal". method : str, optional Method for the calculation, by default "minimize". Options are "minimize" or "least_squares". **kwargs : dict Additional arguments for the calculation. - eos_model: Equation of state model to use for the calculation. Options are "SRK" or "PR". - activity_model: Activity model to use for the calculation. Options are "NRTL" or "UNIFAC". - message: Optional message to display during the calculation.

Returns

dict Equilibrium state of the reaction system.

Notes

Inputs must contain mole_fraction or mole, temperature, and pressure as follows:

• mole_fraction: dict, initial mole fraction of the components in the system.
• mole: dict, initial mole of the components in the system.
• temperature: list, temperature in the form of [value, unit].
• pressure: list, pressure in the form of [value, unit].

For non-ideal solution, the rules of the activity model should be followed.

Activity inputs for NRTL model should contain: - tau: dict, tau values for the components. - alpha: dict, alpha values for the components.

Activity inputs for UNIQUAC model should contain: - q: dict, q values for the components. - r: dict, r values for the components. - tau: dict, tau values for the components.

For both models, dU or a, b, c, and d are required to calculate tau in case tau is not provided. All values should be introduced in the inputs as:

inputs = {
    'mole': mole,
    'temperature': [100, "C"],
    'pressure': [1.0, "bar"],
    'tau': tau,
    'alpha': alpha, ...
    }

Source code in PyReactLab/docs/reactionsystem.py

def equilibrium(
    self,
    inputs: Dict[str, Any],
    conversion: Optional[List[str]] = None,
    gas_mixture: Literal[
        "ideal", "non-ideal"
    ] = "ideal",
    solution: Literal[
        "ideal", "non-ideal"
    ] = "ideal",
    method: Literal[
        'minimize', 'least_squares'
    ] = 'minimize',
    **kwargs
):
    """
    Calculate the equilibrium state of the reaction system.

    Parameters
    ----------
    inputs : dict
        Inputs for the equilibrium calculation.
    conversion : list, optional
        List of components to calculate conversion for, by default None.
    gas_mixture : str, optional
        Type of gas mixture, by default "ideal".
    solution : str, optional
        Type of liquid mixture, by default "ideal".
    method : str, optional
        Method for the calculation, by default "minimize".
        Options are "minimize" or "least_squares".
    **kwargs : dict
        Additional arguments for the calculation.
            - eos_model: Equation of state model to use for the calculation. Options are "SRK" or "PR".
            - activity_model: Activity model to use for the calculation. Options are "NRTL" or "UNIFAC".
            - message: Optional message to display during the calculation.


    Returns
    -------
    dict
        Equilibrium state of the reaction system.

    Notes
    -----
    Inputs must contain `mole_fraction` or `mole`, `temperature`, and `pressure` as follows:

    - mole_fraction: dict, initial mole fraction of the components in the system.
    - mole: dict, initial mole of the components in the system.
    - temperature: list, temperature in the form of [value, unit].
    - pressure: list, pressure in the form of [value, unit].

    For `non-ideal solution`, the rules of the activity model should be followed.

    Activity inputs for NRTL model should contain:
    - tau: dict, tau values for the components.
    - alpha: dict, alpha values for the components.

    Activity inputs for UNIQUAC model should contain:
    - q: dict, q values for the components.
    - r: dict, r values for the components.
    - tau: dict, tau values for the components.

    For both models, dU or a, b, c, and d are required to calculate tau in case tau is not provided. All values should be introduced in the `inputs` as:

    ```python
    inputs = {
        'mole': mole,
        'temperature': [100, "C"],
        'pressure': [1.0, "bar"],
        'tau': tau,
        'alpha': alpha, ...
        }
    ```
    """
    try:
        # ! Start timing
        start_time = time.time()

        # SECTION: check args
        # NOTE: check if inputs are valid
        if not isinstance(inputs, dict):
            raise ValueError("Inputs must be a dictionary.")
        # check if inputs are valid
        if not any(key in inputs for key in ["mole_fraction", "mole"]):
            raise ValueError(
                "Inputs must contain either mole_fraction or mole.")
        if "temperature" not in inputs:
            raise ValueError("Inputs must contain temperature.")
        if "pressure" not in inputs:
            raise ValueError("Inputs must contain pressure.")

        # NOTE: check if temperature is valid
        # set
        temperature = inputs["temperature"]
        # check if temperature is valid
        if not isinstance(temperature, list):
            raise ValueError("Temperature must be a number.")

        # check if temperature is valid
        if len(temperature) != 2:
            raise ValueError(
                "Temperature must be a list of length 2 following [value, unit].")

        if not isinstance(temperature[0], (int, float)):
            raise ValueError("Temperature must be a number.")

        if not isinstance(temperature[1], str):
            raise ValueError("Temperature unit must be a string.")

        # ! convert to **K**
        # set unit
        unit_set = f"{temperature[1]} => K"
        temperature_K = pycuc.to(temperature[0], unit_set)

        # NOTE: check if pressure is valid
        # set
        pressure = inputs["pressure"]
        # check if pressure is valid
        if not isinstance(pressure, list):
            raise ValueError("Pressure must be a number.")

        # check if pressure is valid
        if len(pressure) != 2:
            raise ValueError(
                "Pressure must be a list of length 2 following [value, unit].")

        if not isinstance(pressure[0], (int, float)):
            raise ValueError("Pressure must be a number.")

        if not isinstance(pressure[1], str):
            raise ValueError("Pressure unit must be a string.")

        # ! convert to **bar**
        # set unit
        unit_set = f"{pressure[1]} => bar"
        pressure_bar = pycuc.to(pressure[0], unit_set)

        # NOTE: check if initial mole fraction is valid
        # set
        initial_mole_fraction = inputs.get("mole_fraction", None)
        # check
        initial_mole = inputs.get("mole", None)

        if initial_mole_fraction is not None:
            # check if initial mole fraction is valid
            if not isinstance(initial_mole_fraction, dict):
                raise ValueError(
                    "Initial mole fraction must be a dictionary.")

            # check if initial mole fraction is valid
            for key in initial_mole_fraction:
                if not isinstance(key, str):
                    raise ValueError(
                        "Initial mole fraction key must be a string.")
                if not isinstance(initial_mole_fraction[key], (int, float)):
                    raise ValueError(
                        "Initial mole fraction value must be a number.")

            # ? normalize
            # check if sum of mole fraction is 1
            initial_mole_fraction = ReactionAnalyzer.norm_mole_fraction(
                initial_mole_fraction)

            # NOTE: convert to mole fraction
            initial_mole = ReactionAnalyzer.cal_mole(
                initial_mole_fraction)

        elif initial_mole is not None:
            # check if initial mole is valid
            if not isinstance(initial_mole, dict):
                raise ValueError("Initial mole must be a dictionary.")

            # check if initial mole is valid
            for key in initial_mole:
                if not isinstance(key, str):
                    raise ValueError(
                        "Initial mole key must be a string.")
                if not isinstance(initial_mole[key], (int, float)):
                    raise ValueError(
                        "Initial mole value must be a number.")

            # NOTE: convert to mole fraction
            initial_mole_fraction, _ = ReactionAnalyzer.cal_mole_fraction(
                initial_mole)

        # NOTE: build mole and mole fraction matrix regarding the component dict
        # check
        if (initial_mole is not None and
                initial_mole_fraction is not None):
            # set values
            (
                initial_mole_std, initial_mole_fraction_std, _, _
            ) = ReactionAnalyzer.set_stream(
                component_dict=self.component_dict,
                mole=initial_mole,
                mole_fraction=initial_mole_fraction,
            )
        else:
            raise ValueError(
                "Initial mole and mole fraction must be provided.")

        # SECTION: kwargs
        # eos model (name)
        eos_model = kwargs.get("eos_model", "SRK")
        # check if eos model is valid
        if eos_model not in ["SRK", "PR", 'RK', 'vdW']:
            raise ValueError(
                "Invalid eos model. Options are 'SRK','RK','PR', or 'vdW'.")

        # activity model (name)
        activity_model = kwargs.get("activity_model", "NRTL")
        # check if activity model is valid
        if activity_model not in ["NRTL", "UNIQUAC"]:
            raise ValueError(
                "Invalid activity model. Options are 'NRTL' or 'UNIQUAC'.")

        # SECTION: init
        ReactionOptimizer_ = ReactionOptimizer(
            self.datasource,
            self.equationsource,
            self.component_dict,
            self.coeff_list_dict,
            self.coeff_T_list_list,
            self.reaction_analysis,
            self.overall_reaction_analysis,
        )

        # SECTION: setting up the reaction optimizer
        # NOTE: set up eos model
        # ? eos model
        ReactionOptimizer_.eos_model = eos_model
        # init eos class
        ReactionOptimizer_.eos = ptm.eos()

        # gas mixture
        ReactionOptimizer_.gas_mixture = gas_mixture

        # NOTE: set up activity model
        # ? activity inputs
        # NRTL: tau, alpha,
        # dg or a, b, c, and d are required to calculated tau
        # UNIQUAC: q, r, tau
        # dU or a, b, c, and d are required to calculated tau
        activity_inputs = {**inputs}
        # remove temperature and pressure
        activity_inputs.pop("temperature", None)
        activity_inputs.pop("pressure", None)
        # remove mole and mole fraction
        activity_inputs.pop("mole", None)
        activity_inputs.pop("mole_fraction", None)

        # set
        if activity_inputs is not None:
            #  check if activity inputs is valid
            if not isinstance(activity_inputs, dict):
                raise ValueError(
                    "Activity inputs must be a dictionary.")

            # set
            ReactionOptimizer_.activity_inputs = activity_inputs

        # ? activity model
        ReactionOptimizer_.activity_model = activity_model
        # init activity class
        ReactionOptimizer_.activity = ptm.activities(
            components=self.component_list,
            model_name=activity_model,
            model_source=self.model_source,
        )

        # solution
        ReactionOptimizer_.solution = solution

        # SECTION: equilibrium constant calculation
        # res
        equilibrium_constant = {}

        # loop through each reaction
        for key, r in self.__reaction_list.items():
            # NOTE: check if reaction is valid
            if not isinstance(r, Reaction):
                raise ValueError(
                    f"Invalid reaction object for {key}")

            # ! calculate equilibrium constant at the given temperature
            res_ = r.cal_equilibrium_constant(temperature)

            # save
            equilibrium_constant[key] = res_

        # NOTE: run equilibrium calculation
        opt_res = ReactionOptimizer_.opt_run(
            initial_mole=initial_mole_std,
            initial_mole_fraction=initial_mole_fraction_std,
            temperature=temperature_K,
            pressure=pressure_bar,
            equilibrium_constants=equilibrium_constant,
            reaction_numbers=self.reaction_numbers,
            method=method,
        )

        # NOTE: check optimization result
        res = {}

        # check if optimization is successful
        if opt_res.success:
            # extent of reaction
            EoR = opt_res.x

            # calculate equilibrium state
            Eq_Xfs, Eq_Xfs_vector, Eq_Nf, Eq_Nfs_vector, Eq_Nfs = \
                ReactionOptimizer_.equilibrium_results(
                    initial_mole=initial_mole_std,
                    EoR=EoR,
                )

            # NOTE: calculate conversion
            # res
            conversion_res = None
            # conversion
            if conversion is not None:
                # check list
                if not isinstance(conversion, list):
                    raise ValueError(
                        "Conversion must be a list of strings.")

                # check component in the components
                for key in conversion:
                    if key not in self.component_dict:
                        raise ValueError(
                            f"Invalid component: {key} in the components.")

                # calculate conversion
                conversion_res = ReactionAnalyzer.cal_conversion(
                    initial_mole=initial_mole_std,
                    final_mole=Eq_Nfs,
                    components=conversion,
                )

            # SECTION: set results
            # initial feed
            res['feed'] = {
                "mole": {
                    'value': initial_mole_std,
                    'unit': "mol"
                },
                "mole_total": {
                    'value': sum(initial_mole_std.values()),
                    'unit': "mol"
                },
                "mole_fraction": {
                    'value': initial_mole_fraction_std,
                    'unit': "dimensionless"
                },
                "mole_fraction_sum": {
                    'value': sum(initial_mole_fraction_std.values()),
                    'unit': "dimensionless"
                },
            }

            # equilibrium state
            res['equilibrium'] = {
                "mole": {
                    'value': Eq_Nfs,
                    'unit': "mol"
                },
                "mole_total": {
                    'value': sum(Eq_Nfs.values()),
                    'unit': "mol"
                },
                "mole_fraction": {
                    'value': Eq_Xfs,
                    'unit': "dimensionless"
                },
                "mole_fraction_sum": {
                    'value': sum(Eq_Xfs.values()),
                    'unit': "dimensionless"
                },
            }

            # extent of reaction
            res['extent_of_reaction'] = {
                "value": EoR,
                "unit": "dimensionless"
            }

            # equilibrium condition
            # temperature
            res['temperature'] = {
                "value": temperature_K,
                "unit": "K"
            }

            # pressure
            res['pressure'] = {
                "value": pressure_bar,
                "unit": "bar"
            }

            # optimization fun value
            res['optimization_fun'] = {
                "value": opt_res.fun,
                "unit": "dimensionless"
            }

            # conversion
            if conversion_res is not None:
                res['conversion'] = conversion_res

            # message
            message = kwargs.get("message", None)
            if message is not None:
                res['message'] = message

            # gas mixture
            res['gas_mixture'] = gas_mixture
            # solution
            res['solution'] = solution

        # NOTE: set time
        # ! Stop timing
        end_time = time.time()
        computation_time = end_time - start_time
        # add to res
        res['computation_time'] = {
            "value": computation_time,
            "unit": "s"
        }

        # res
        return res

    except Exception as e:
        raise Exception(
            f"Failing in the equilibrium calculations {str(e)}") from e

equilibrium_deviation(inputs, gas_mixture='ideal', solution='ideal', **kwargs)

Conduct equilibrium deviation calculations for the reaction system.

This method calculates the actual Gibbs free energy of reaction (GiEn_rxn_T) at a given temperature and pressure, and component composition.

There are three cases as: - Case 1: GiEn_rxn_T < 0, the reaction is not at equilibrium and the reaction will proceed in the forward direction. - Case 2: GiEn_rxn_T > 0, the reaction is not at equilibrium and the reaction will proceed in the reverse direction. - Case 3: GiEn_rxn_T = 0, the reaction is at equilibrium and no further reaction will occur.

Parameters

inputs : dict Inputs for the equilibrium deviation calculation. conversion : list, optional List of components to calculate conversion for, by default None. gas_mixture : str, optional Type of gas mixture, by default "ideal". solution : str, optional Type of liquid mixture, by default "ideal". **kwargs : dict Additional arguments for the calculation. - eos_model: Equation of state model to use for the calculation. Options are "SRK" or "PR". - activity_model: Activity model to use for the calculation. Options are "NRTL" or "UNIFAC". - message: Optional message to display during the calculation. - mole_basis: 1 mole - minimum_mole: Minimum mole for the initial mole, by default 1e-5.

Returns

dict Equilibrium deviation of the reaction system.

Notes

The inputs and kwargs are similar to those in the equilibrium method.

Source code in PyReactLab/docs/reactionsystem.py

def equilibrium_deviation(
    self,
    inputs: Dict[str, Any],
    gas_mixture: Literal[
        "ideal", "non-ideal"
    ] = "ideal",
    solution: Literal[
        "ideal", "non-ideal"
    ] = "ideal",
    **kwargs
):
    """
    Conduct equilibrium deviation calculations for the reaction system.

    This method calculates `the actual Gibbs free energy of reaction` (GiEn_rxn_T) at a given temperature and pressure, and component composition.

    There are three cases as:
    - Case 1: `GiEn_rxn_T < 0`, the reaction is not at equilibrium and the reaction will proceed in the forward direction.
    - Case 2: `GiEn_rxn_T > 0`, the reaction is not at equilibrium and the reaction will proceed in the reverse direction.
    - Case 3: `GiEn_rxn_T = 0`, the reaction is at equilibrium and no further reaction will occur.

    Parameters
    ----------
    inputs : dict
        Inputs for the equilibrium deviation calculation.
    conversion : list, optional
        List of components to calculate conversion for, by default None.
    gas_mixture : str, optional
        Type of gas mixture, by default "ideal".
    solution : str, optional
        Type of liquid mixture, by default "ideal".
    **kwargs : dict
        Additional arguments for the calculation.
            - eos_model: Equation of state model to use for the calculation. Options are "SRK" or "PR".
            - activity_model: Activity model to use for the calculation. Options are "NRTL" or "UNIFAC".
            - message: Optional message to display during the calculation.
            - mole_basis: 1 mole
            - minimum_mole: Minimum mole for the initial mole, by default 1e-5.

    Returns
    -------
    dict
        Equilibrium deviation of the reaction system.

    Notes
    -----
    The inputs and kwargs are similar to those in the `equilibrium` method.
    """
    # NOTE: call equilibrium method with deviation=True
    try:
        # ! Start timing
        start_time = time.time()

        if not self.overall_reaction_phase:
            raise ValueError(
                "Overall reaction phase is not set. Please run the primary analysis first.")

        # mole basis
        mole_basis = kwargs.get("mole_basis", 1.0)
        # minimum mole
        minimum_mole = kwargs.get("minimum_mole", 1e-5)

        # SECTION: check args
        # NOTE: check if inputs are valid
        if not isinstance(inputs, dict):
            raise ValueError("Inputs must be a dictionary.")
        # check if inputs are valid
        if not any(key in inputs for key in ["mole_fraction", "mole"]):
            raise ValueError(
                "Inputs must contain either mole_fraction or mole.")
        if "temperature" not in inputs:
            raise ValueError("Inputs must contain temperature.")
        if "pressure" not in inputs:
            raise ValueError("Inputs must contain pressure.")

        # SECTION: set temperature and pressure
        # NOTE: check if temperature is valid
        # set
        temperature = inputs["temperature"]
        # check if temperature is valid
        if not isinstance(temperature, list):
            raise ValueError("Temperature must be a number.")

        # check if temperature is valid
        if len(temperature) != 2:
            raise ValueError(
                "Temperature must be a list of length 2 following [value, unit].")

        if not isinstance(temperature[0], (int, float)):
            raise ValueError("Temperature must be a number.")

        if not isinstance(temperature[1], str):
            raise ValueError("Temperature unit must be a string.")

        # ! convert to **K**
        # set unit
        unit_set = f"{temperature[1]} => K"
        temperature_K = pycuc.to(temperature[0], unit_set)
        # set temperature
        temperature = Temperature(
            value=temperature_K,
            unit="K"
        )

        # NOTE: check if pressure is valid
        # set
        pressure = inputs["pressure"]
        # check if pressure is valid
        if not isinstance(pressure, list):
            raise ValueError("Pressure must be a number.")

        # check if pressure is valid
        if len(pressure) != 2:
            raise ValueError(
                "Pressure must be a list of length 2 following [value, unit].")

        if not isinstance(pressure[0], (int, float)):
            raise ValueError("Pressure must be a number.")

        if not isinstance(pressure[1], str):
            raise ValueError("Pressure unit must be a string.")

        # ! convert to **bar**
        # set unit
        unit_set = f"{pressure[1]} => bar"
        pressure_bar = pycuc.to(pressure[0], unit_set)
        # set pressure
        pressure = Pressure(
            value=pressure_bar,
            unit="bar"
        )

        # NOTE: set operating conditions
        operating_conditions = OperatingConditions(
            temperature=temperature,
            pressure=pressure,
        )

        # SECTION:
        # NOTE: check if initial mole fraction is valid
        # set
        initial_mole_fraction = inputs.get("mole_fraction", None)
        # check
        initial_mole = inputs.get("mole", None)

        if initial_mole_fraction is not None:
            # check if initial mole fraction is valid
            if not isinstance(initial_mole_fraction, dict):
                raise ValueError(
                    "Initial mole fraction must be a dictionary.")

            # check if initial mole fraction is valid
            for key in initial_mole_fraction:
                if not isinstance(key, str):
                    raise ValueError(
                        "Initial mole fraction key must be a string.")
                if not isinstance(initial_mole_fraction[key], (int, float)):
                    raise ValueError(
                        "Initial mole fraction value must be a number.")

            # check if sum of mole fraction is 1
            initial_mole_fraction = ReactionAnalyzer.norm_mole_fraction(
                initial_mole_fraction)

            # NOTE: convert to mole fraction
            initial_mole = ReactionAnalyzer.cal_mole(
                initial_mole_fraction,
                mole_basis
            )

        elif initial_mole is not None:
            # check if initial mole is valid
            if not isinstance(initial_mole, dict):
                raise ValueError("Initial mole must be a dictionary.")

            # check if initial mole is valid
            for key in initial_mole:
                if not isinstance(key, str):
                    raise ValueError(
                        "Initial mole key must be a string.")
                if not isinstance(initial_mole[key], (int, float)):
                    raise ValueError(
                        "Initial mole value must be a number.")

            # set initial mole
            initial_mole = ReactionAnalyzer.set_initial_mole(
                initial_mole,
                minimum_mole=minimum_mole
            )

            # NOTE: convert to mole fraction
            initial_mole_fraction, _ = ReactionAnalyzer.cal_mole_fraction(
                initial_mole)

        # NOTE: build mole and mole fraction matrix regarding the component dict
        # check
        if initial_mole is not None and initial_mole_fraction is not None:
            # set values
            (
                initial_mole_std, initial_mole_fraction_std, _, _
            ) = ReactionAnalyzer.set_stream(
                component_dict=self.component_dict,
                mole=initial_mole,
                mole_fraction=initial_mole_fraction,
            )
        else:
            raise ValueError(
                "Initial mole and mole fraction must be provided.")

        # SECTION: mole and mole fraction for gas mixture and solution
        self.phase_stream = ReactionAnalyzer.set_phase_stream(
            initial_mole=initial_mole_std,
            initial_mole_fraction=initial_mole_fraction_std,
            phase_contents=self.phase_contents,
        )

        # SECTION: kwargs
        # eos model (name)
        eos_model = kwargs.get("eos_model", "SRK")
        # check if eos model is valid
        if eos_model not in ["SRK", "PR", 'RK', 'vdW']:
            raise ValueError(
                "Invalid eos model. Options are 'SRK','RK','PR', or 'vdW'.")

        # activity model (name)
        activity_model = kwargs.get("activity_model", "NRTL")
        # check if activity model is valid
        if activity_model not in ["NRTL", "UNIQUAC"]:
            raise ValueError(
                "Invalid activity model. Options are 'NRTL' or 'UNIQUAC'.")

        # SECTION: init chemical potential
        # init
        ChemicalPotential_ = ChemicalPotential(
            self.datasource,
            self.equationsource,
            self.reaction_list,
            self.component_dict,
            self.component_state_list,
            self.coeff_list_dict,
            self.reaction_analysis,
            self.phase_stream,
            self.phase_contents,
            self.overall_reaction_analysis,
            self.overall_reaction_phase,
            operating_conditions=operating_conditions,
        )

        # SECTION: setting up the chemical potential
        # NOTE: set up eos model
        # ? eos model
        ChemicalPotential_.eos_model = eos_model
        # init eos class
        ChemicalPotential_.eos = ptm.eos()

        # gas mixture
        ChemicalPotential_.gas_mixture = gas_mixture

        # NOTE: set up activity model
        # ? activity inputs
        # NRTL: tau, alpha,
        # dg or a, b, c, and d are required to calculated tau
        # UNIQUAC: q, r, tau
        # dU or a, b, c, and d are required to calculated tau
        activity_inputs = {**inputs}
        # remove temperature and pressure
        activity_inputs.pop("temperature", None)
        activity_inputs.pop("pressure", None)
        # remove mole and mole fraction
        activity_inputs.pop("mole", None)
        activity_inputs.pop("mole_fraction", None)

        # set
        if activity_inputs is not None:
            #  check if activity inputs is valid
            if not isinstance(activity_inputs, dict):
                raise ValueError(
                    "Activity inputs must be a dictionary.")

            # set
            ChemicalPotential_.activity_inputs = activity_inputs

        # ? activity model
        ChemicalPotential_.activity_model = activity_model
        # ? set component list (liquid phase)
        # components_l = self.phase_contents['l']

        # init activity class
        ChemicalPotential_.activity = ptm.activities(
            components=self.component_list,
            model_name=activity_model,
            model_source=self.model_source,
        )

        # solution
        ChemicalPotential_.solution = solution

        # SECTION: the actual Gibbs free energy of reaction
        # calc
        res = ChemicalPotential_.cal_actual_gibbs_energy_of_reaction(
            operating_conditions=operating_conditions,
            gas_mixture=gas_mixture,
            solution=solution,
        )

        # NOTE: set time
        # ! Stop timing
        end_time = time.time()
        computation_time = end_time - start_time
        # add to res
        res['computation_time'] = {
            "value": computation_time,
            "unit": "s"
        }

        # res
        return res
    except Exception as e:
        raise Exception(
            f"Failing in the equilibrium deviation calculations {str(e)}") from e

reaction_equilibrium_constant(reaction_name, temperature, method="van't Hoff", **kwargs)

Calculate the equilibrium constant at a given temperature using the van't Hoff equation.

Parameters

reaction_name : str Name of the reaction. temperature : list Temperature in the form of [value, unit], the unit is automatically converted to K. method : str, optional Method to calculate the equilibrium constant, by default "van't Hoff". Options are "van't Hoff" or "shortcut van't Hoff". **kwargs : dict Additional arguments for the calculation. - message: Optional message to display during the calculation.

Returns

res : dict Equilibrium constant at the given temperature.

Source code in PyReactLab/docs/reactionsystem.py

def reaction_equilibrium_constant(
    self,
    reaction_name: str,
    temperature: list[float | str],
    method: Literal[
        "van't Hoff", "shortcut van't Hoff"
    ] = "van't Hoff",
    **kwargs
):
    """
    Calculate the equilibrium constant at a given temperature using the van't Hoff equation.

    Parameters
    ----------
    reaction_name : str
        Name of the reaction.
    temperature : list
        Temperature in the form of [value, unit], the unit is automatically converted to K.
    method : str, optional
        Method to calculate the equilibrium constant, by default "van't Hoff".
        Options are "van't Hoff" or "shortcut van't Hoff".
    **kwargs : dict
        Additional arguments for the calculation.
            - message: Optional message to display during the calculation.

    Returns
    -------
    res : dict
        Equilibrium constant at the given temperature.
    """
    try:
        # NOTE: check if temperature is valid
        if not isinstance(temperature, list):
            raise ValueError("Temperature must be a number.")

        # check if temperature is valid
        if len(temperature) != 2:
            raise ValueError(
                "Temperature must be a list of length 2 following [value, unit].")

        if not isinstance(temperature[0], (int, float)):
            raise ValueError("Temperature must be a number.")

        if not isinstance(temperature[1], str):
            raise ValueError("Temperature unit must be a string.")

        # NOTE: check if reaction name is valid
        if reaction_name not in self.__reaction_list:
            raise ValueError(f"Invalid reaction name: {reaction_name}")

        # SECTION: pre-calculation
        # get reaction object
        reaction = self.__reaction_list[reaction_name]

        # check if reaction is valid
        if not isinstance(reaction, Reaction):
            raise ValueError(
                f"Invalid reaction object for {reaction_name}")

        # NOTE: calculate equilibrium constant at the given temperature
        res = reaction.cal_equilibrium_constant(
            temperature,
            method=method
        )

        # NOTE: message
        message = kwargs.get("message", None)
        if message is not None:
            res['message'] = message

        # res
        return res
    except Exception as e:
        raise Exception(
            f"Error in ReactionSystem.equilibrium_constant_at_temperature(): {str(e)}") from e

select_reaction(reaction_name)

Select a reaction from the reaction system.

Parameters

reaction_name : str Name of the reaction.

Returns

Reaction or None Reaction object if found, otherwise None.

Source code in PyReactLab/docs/reactionsystem.py

def select_reaction(self, reaction_name: str) -> Reaction:
    """
    Select a reaction from the reaction system.

    Parameters
    ----------
    reaction_name : str
        Name of the reaction.

    Returns
    -------
    Reaction or None
        Reaction object if found, otherwise None.
    """
    try:
        # NOTE: check if reaction name is valid
        if reaction_name not in self.__reaction_list:
            raise ValueError(f"Invalid reaction name: {reaction_name}")

        # reaction object
        reaction = self.__reaction_list.get(reaction_name, None)
        # check if reaction is valid
        if reaction is None or reaction == 'None':
            raise ValueError(
                f"Invalid reaction object for {reaction_name}")
        return reaction
    except Exception as e:
        raise Exception(
            f"Error in ReactionSystem.select_reaction(): {str(e)}") from e

🔬 Reaction

Reaction

Class to represent a chemical reaction.

Source code in PyReactLab/docs/reaction.py

class Reaction():
    """Class to represent a chemical reaction."""
    # NOTE: variables

    # variables
    reaction_analysis_result = {}

    def __init__(self,
                 datasource,
                 equationsource,
                 reaction,
                 **kwargs):
        '''
        Initialize the reaction object.

        Parameters
        ----------
        datasource : dict
            Dictionary containing the thermodynamic data for the reaction.
        equationsource : dict
            Dictionary containing the equations for the reaction.
        reaction : dict
            Reaction object in the form of a dictionary with the following keys:
            - 'name': str, the name of the reaction.
            - 'reaction': str, the reaction equation.
        **kwargs : dict
            Additional keyword arguments.
        '''
        # set
        self.datasource = datasource
        self.equationsource = equationsource
        self.reaction = reaction

        # NOTE: init
        self.ChemReactUtils_ = ChemReactUtils()
        self.ReactionAnalyzer_ = ReactionAnalyzer()

        # NOTE: kwargs
        # phase rule
        self.phase_rule = kwargs.get('phase_rule', None)

        # NOTE: reaction analyzer
        self._reaction_analyzer()

    def _reaction_analyzer(self) -> None:
        """
        Execute the primary analysis for the reaction system.
        """
        try:
            # SECTION: reaction system analysis
            # analyze reaction
            reaction_res = {}

            # NOTE: name
            name = self.reaction['name']
            # update
            reaction_res['name'] = name

            # NOTE: reaction
            reaction = self.reaction['reaction']
            # update
            reaction_res['reaction'] = reaction

            # NOTE: extract reaction information
            # 'name': name,
            # 'reaction': reaction,
            # 'reactants': reactants,
            # 'products': products,
            # 'reaction_coefficient': reaction_coefficient,
            # 'carbon_count': carbon_count
            # 'reaction_state': reaction_state,
            _res_0 = self.ChemReactUtils_.analyze_reaction(
                self.reaction,
                phase_rule=self.phase_rule,
            )

            # NOTE: energy analysis
            # set input
            _input = {**reaction_res, **_res_0}

            energy_analysis = self.ReactionAnalyzer_.energy_analysis(
                self.datasource,
                self.equationsource,
                _input
            )

            # update
            self.reaction_name = name
            self.reaction_equation = reaction
            self.reactants = _res_0['reactants']
            self.reactants_names = _res_0['reactants_names']
            self.products = _res_0['products']
            self.products_names = _res_0['products_names']
            self.reaction_coefficient = _res_0['reaction_coefficient']
            self.carbon_count = _res_0['carbon_count']
            self.reaction_state = _res_0['reaction_state']
            self.reaction_phase = _res_0['reaction_phase']
            self.state_count = _res_0['state_count']
            self.energy_analysis = energy_analysis
            # extract
            self.gibbs_free_energy_of_formation_std = energy_analysis[
                GIBBS_FREE_ENERGY_OF_FORMATION_STD]
            self.enthalpy_of_formation_std = energy_analysis[
                ENTHALPY_OF_FORMATION_STD
            ]
            self.gibbs_free_energy_of_reaction_std = energy_analysis[
                GIBBS_FREE_ENERGY_OF_REACTION_STD]
            self.enthalpy_of_reaction_std = energy_analysis[
                ENTHALPY_OF_REACTION_STD
            ]
            self.equilibrium_constant_std = energy_analysis[
                EQUILIBRIUM_CONSTANT_STD
            ]

            # reaction analysis result
            self.reaction_analysis_result: Dict[str, Any] = {
                'name': name,
                'reaction': reaction,
                'reactants': _res_0['reactants'],
                'reactants_names': _res_0['reactants_names'],
                'products': _res_0['products'],
                'products_names': _res_0['products_names'],
                'reaction_coefficient': _res_0['reaction_coefficient'],
                'carbon_count': _res_0['carbon_count'],
                'reaction_state': _res_0['reaction_state'],
                'state_count': _res_0['state_count'],
                'energy_analysis': energy_analysis
            }

        except Exception as e:
            raise Exception(
                f"Error in ReactionSystem.go(): {str(e)}") from e

    def check_state(self):
        """
        Check the state of the reaction system whether it is a liquid, gas, aqueous, or solid or a combination of them.
        """
        try:
            # NOTE: set state
            states = []

            # looping through reactants and products
            for reactant in self.reactants:
                _res = reactant['state']
                # ? add state
                states.append(_res)
            for product in self.products:
                _res = product['state']
                # ? add state
                states.append(_res)

            # NOTE: remove duplicates
            states = list(set(states))
        except Exception as e:
            raise Exception(
                f"Failing in finding the reaction state: {str(e)}") from e

    def cal_equilibrium_constant(
        self,
        temperature: List[float | str],
        method: Literal[
            "van't Hoff", "shortcut van't Hoff"
        ] = "van't Hoff"
    ) -> Dict[str, Any]:
        """
        Calculate the equilibrium constant at a given temperature using the van't Hoff equation.

        Parameters
        ----------
        temperature : list[float, str]
            Temperature in any unit, e.g. [300.0, "K"], It is automatically converted to Kelvin.
        method : str, optional
            Method to calculate the equilibrium constant. The default is "van't Hoff".
            Options are "van't Hoff" or "shortcut van't Hoff".

        Returns
        -------
        dict
            Dictionary containing the equilibrium constant at the given temperature.
        """
        try:
            # NOTE: check if T
            # check if T is a list
            if not isinstance(temperature, list):
                raise ValueError("Temperature must be a list.")

            # check if T is a number
            if not isinstance(temperature[0], (int, float)):
                raise ValueError("Temperature must be a number.")

            # check if T is a string
            if not isinstance(temperature[1], str):
                raise ValueError("Temperature unit must be a string.")

            # NOTE: convert temperature to Kelvin
            # set unit
            unit_set = f"{temperature[1]} => K"
            T = pycuc.to(temperature[0], unit_set)

            # SECTION: calculate equilibrium constant at a given temperature
            # check if method is valid
            if method not in ["van't Hoff", "shortcut van't Hoff"]:
                raise ValueError(
                    f"Invalid method: {method}, must be 'van't Hoff' or 'shortcut van't Hoff'.")

            if method == "van't Hoff":
                # ! NOTE: van't Hoff method
                res = self.ReactionAnalyzer_.vh(
                    self.datasource,
                    self.equationsource,
                    T,
                    self.reaction_analysis_result,
                )

                # update
                res['method'] = method
                # return
                return res
            elif method == "shortcut van't Hoff":
                # ! NOTE: shortcut van't Hoff method
                res = self.ReactionAnalyzer_.vh_shortcut(
                    self.datasource,
                    self.equationsource,
                    T,
                    self.enthalpy_of_reaction_std['value'],
                    self.equilibrium_constant_std['value'],
                    self.reaction_name,
                    self.reaction_equation
                )

                # update
                res['method'] = method
                # return
                return res
            else:
                raise ValueError(
                    f"Invalid method: {method}, must be 'van't Hoff' or 'shortcut van't Hoff'.")

        except Exception as e:
            raise Exception(
                f"Error in ReactionSystem.Keq_T(): {str(e)}") from e

    def cal_reaction_energy(
        self,
        temperature: List[float | str]
    ) -> Dict[str, Any]:
        """
        Calculate the reaction energy at a given temperature which consists of the following:
        - Gibbs free energy of reaction at T
        - Enthalpy of reaction at T

        Parameters
        ----------
        temperature : list[float, str]
            Temperature in any unit, e.g. [300.0, "K"], It is automatically converted to Kelvin.

        Returns
        -------
        dict
            Dictionary containing the reaction energy at the given temperature.
        """
        try:
            # NOTE: check if T
            # check if T is a list
            if not isinstance(temperature, list):
                raise ValueError("Temperature must be a list.")

            # check if T is a number
            if not isinstance(temperature[0], (int, float)):
                raise ValueError("Temperature must be a number.")

            # check if T is a string
            if not isinstance(temperature[1], str):
                raise ValueError("Temperature unit must be a string.")

            # NOTE: convert temperature to Kelvin
            # set unit
            unit_set = f"{temperature[1]} => K"
            T = pycuc.to(temperature[0], unit_set)

            # SECTION: calculate reaction energy at a given temperature
            res = self.ReactionAnalyzer_.reaction_energy_analysis(
                self.datasource,
                self.equationsource,
                T,
                self.reaction_analysis_result,
            )

            # NOTE: extract reaction energy
            return {
                "reaction_name": self.reaction_name,
                "reaction_equation": self.reaction_equation,
                "temperature": {
                    "value": T,
                    "symbol": "T",
                    "unit": "K",
                },
                GIBBS_FREE_ENERGY_OF_REACTION_T: res[GIBBS_FREE_ENERGY_OF_REACTION_T],
                ENTHALPY_OF_REACTION_T: res[ENTHALPY_OF_REACTION_T],
            }

        except Exception as e:
            raise Exception(
                f"Error in ReactionSystem.Keq_T(): {str(e)}") from e

__init__(datasource, equationsource, reaction, **kwargs)

Initialize the reaction object.

Parameters

datasource : dict Dictionary containing the thermodynamic data for the reaction. equationsource : dict Dictionary containing the equations for the reaction. reaction : dict Reaction object in the form of a dictionary with the following keys: - 'name': str, the name of the reaction. - 'reaction': str, the reaction equation. **kwargs : dict Additional keyword arguments.

Source code in PyReactLab/docs/reaction.py

def __init__(self,
             datasource,
             equationsource,
             reaction,
             **kwargs):
    '''
    Initialize the reaction object.

    Parameters
    ----------
    datasource : dict
        Dictionary containing the thermodynamic data for the reaction.
    equationsource : dict
        Dictionary containing the equations for the reaction.
    reaction : dict
        Reaction object in the form of a dictionary with the following keys:
        - 'name': str, the name of the reaction.
        - 'reaction': str, the reaction equation.
    **kwargs : dict
        Additional keyword arguments.
    '''
    # set
    self.datasource = datasource
    self.equationsource = equationsource
    self.reaction = reaction

    # NOTE: init
    self.ChemReactUtils_ = ChemReactUtils()
    self.ReactionAnalyzer_ = ReactionAnalyzer()

    # NOTE: kwargs
    # phase rule
    self.phase_rule = kwargs.get('phase_rule', None)

    # NOTE: reaction analyzer
    self._reaction_analyzer()

cal_equilibrium_constant(temperature, method="van't Hoff")

Calculate the equilibrium constant at a given temperature using the van't Hoff equation.

Parameters

temperature : list[float, str] Temperature in any unit, e.g. [300.0, "K"], It is automatically converted to Kelvin. method : str, optional Method to calculate the equilibrium constant. The default is "van't Hoff". Options are "van't Hoff" or "shortcut van't Hoff".

Returns

dict Dictionary containing the equilibrium constant at the given temperature.

Source code in PyReactLab/docs/reaction.py

def cal_equilibrium_constant(
    self,
    temperature: List[float | str],
    method: Literal[
        "van't Hoff", "shortcut van't Hoff"
    ] = "van't Hoff"
) -> Dict[str, Any]:
    """
    Calculate the equilibrium constant at a given temperature using the van't Hoff equation.

    Parameters
    ----------
    temperature : list[float, str]
        Temperature in any unit, e.g. [300.0, "K"], It is automatically converted to Kelvin.
    method : str, optional
        Method to calculate the equilibrium constant. The default is "van't Hoff".
        Options are "van't Hoff" or "shortcut van't Hoff".

    Returns
    -------
    dict
        Dictionary containing the equilibrium constant at the given temperature.
    """
    try:
        # NOTE: check if T
        # check if T is a list
        if not isinstance(temperature, list):
            raise ValueError("Temperature must be a list.")

        # check if T is a number
        if not isinstance(temperature[0], (int, float)):
            raise ValueError("Temperature must be a number.")

        # check if T is a string
        if not isinstance(temperature[1], str):
            raise ValueError("Temperature unit must be a string.")

        # NOTE: convert temperature to Kelvin
        # set unit
        unit_set = f"{temperature[1]} => K"
        T = pycuc.to(temperature[0], unit_set)

        # SECTION: calculate equilibrium constant at a given temperature
        # check if method is valid
        if method not in ["van't Hoff", "shortcut van't Hoff"]:
            raise ValueError(
                f"Invalid method: {method}, must be 'van't Hoff' or 'shortcut van't Hoff'.")

        if method == "van't Hoff":
            # ! NOTE: van't Hoff method
            res = self.ReactionAnalyzer_.vh(
                self.datasource,
                self.equationsource,
                T,
                self.reaction_analysis_result,
            )

            # update
            res['method'] = method
            # return
            return res
        elif method == "shortcut van't Hoff":
            # ! NOTE: shortcut van't Hoff method
            res = self.ReactionAnalyzer_.vh_shortcut(
                self.datasource,
                self.equationsource,
                T,
                self.enthalpy_of_reaction_std['value'],
                self.equilibrium_constant_std['value'],
                self.reaction_name,
                self.reaction_equation
            )

            # update
            res['method'] = method
            # return
            return res
        else:
            raise ValueError(
                f"Invalid method: {method}, must be 'van't Hoff' or 'shortcut van't Hoff'.")

    except Exception as e:
        raise Exception(
            f"Error in ReactionSystem.Keq_T(): {str(e)}") from e

cal_reaction_energy(temperature)

Calculate the reaction energy at a given temperature which consists of the following: - Gibbs free energy of reaction at T - Enthalpy of reaction at T

Parameters

temperature : list[float, str] Temperature in any unit, e.g. [300.0, "K"], It is automatically converted to Kelvin.

Returns

dict Dictionary containing the reaction energy at the given temperature.

Source code in PyReactLab/docs/reaction.py

def cal_reaction_energy(
    self,
    temperature: List[float | str]
) -> Dict[str, Any]:
    """
    Calculate the reaction energy at a given temperature which consists of the following:
    - Gibbs free energy of reaction at T
    - Enthalpy of reaction at T

    Parameters
    ----------
    temperature : list[float, str]
        Temperature in any unit, e.g. [300.0, "K"], It is automatically converted to Kelvin.

    Returns
    -------
    dict
        Dictionary containing the reaction energy at the given temperature.
    """
    try:
        # NOTE: check if T
        # check if T is a list
        if not isinstance(temperature, list):
            raise ValueError("Temperature must be a list.")

        # check if T is a number
        if not isinstance(temperature[0], (int, float)):
            raise ValueError("Temperature must be a number.")

        # check if T is a string
        if not isinstance(temperature[1], str):
            raise ValueError("Temperature unit must be a string.")

        # NOTE: convert temperature to Kelvin
        # set unit
        unit_set = f"{temperature[1]} => K"
        T = pycuc.to(temperature[0], unit_set)

        # SECTION: calculate reaction energy at a given temperature
        res = self.ReactionAnalyzer_.reaction_energy_analysis(
            self.datasource,
            self.equationsource,
            T,
            self.reaction_analysis_result,
        )

        # NOTE: extract reaction energy
        return {
            "reaction_name": self.reaction_name,
            "reaction_equation": self.reaction_equation,
            "temperature": {
                "value": T,
                "symbol": "T",
                "unit": "K",
            },
            GIBBS_FREE_ENERGY_OF_REACTION_T: res[GIBBS_FREE_ENERGY_OF_REACTION_T],
            ENTHALPY_OF_REACTION_T: res[ENTHALPY_OF_REACTION_T],
        }

    except Exception as e:
        raise Exception(
            f"Error in ReactionSystem.Keq_T(): {str(e)}") from e

check_state()

Check the state of the reaction system whether it is a liquid, gas, aqueous, or solid or a combination of them.

Source code in PyReactLab/docs/reaction.py

def check_state(self):
    """
    Check the state of the reaction system whether it is a liquid, gas, aqueous, or solid or a combination of them.
    """
    try:
        # NOTE: set state
        states = []

        # looping through reactants and products
        for reactant in self.reactants:
            _res = reactant['state']
            # ? add state
            states.append(_res)
        for product in self.products:
            _res = product['state']
            # ? add state
            states.append(_res)

        # NOTE: remove duplicates
        states = list(set(states))
    except Exception as e:
        raise Exception(
            f"Failing in finding the reaction state: {str(e)}") from e

🛠️ Optimizer

ReactionOptimizer

Reaction Optimizer class

Source code in PyReactLab/docs/optim.py

class ReactionOptimizer:
    '''Reaction Optimizer class'''

    # SECTION: variables
    _system_inputs = None
    # universal gas constant [J/mol.K]
    R = R_CONST_J__molK
    # temperature [K]
    T_Ref_K = TEMPERATURE_REF_K
    # pressure [Pa]
    P_Ref_Pa = PRESSURE_REF_Pa
    # pressure [bar]
    P_Ref_bar = 1

    # NOTE: models
    # eos
    _eos = None
    # eos model
    _eos_model = None
    # activity model
    _activity_model = None
    # activity
    _activity = None
    # activity inputs
    _activity_inputs = None

    # NOTE: systems
    # gas phase
    _gas_mixture = 'ideal'
    # solution
    _solution = 'ideal'

    def __init__(
        self,
        datasource: Dict[str, Any],
        equationsource: Dict[str, Any],
        component_dict: Dict[str, float | int],
        comp_list: List[Dict[str, float | int]],
        stoichiometric_coeff: np.ndarray,
        reaction_analysis: Dict,
        overall_reaction_analysis: Dict,
        **kwargs
    ):
        '''
        Initialize ReactionOptimizer class

        Parameters
        ----------
        datasource : dict
            datasource dictionary
        equationsource : dict
            equationsource dictionary
        component_dict : dict
            component dictionary {key: value}, such as {CO2: 0, H2: 1, CO: 2, H2O: 3, CH3OH: 4}
        comp_list : list
            component list [{key: value}, {key: value}, ...], such as [{CO2: -1.0, H2: -3.0, CH3OH: 1.0, H2O: 1.0}, ...]
        reaction_analysis : dict
            reaction analysis result
        overall_reaction_analysis : dict
            overall reaction analysis result
        kwargs : dict
            additional parameters
        '''
        # datasource
        self.datasource = datasource
        # equationsource
        self.equationsource = equationsource
        # set component dictionary
        self.component_dict = component_dict
        # set component list
        self.comp_list = comp_list
        # set stoichiometric coefficient
        self.stoichiometric_coeff = stoichiometric_coeff
        # set reaction analysis
        self.reaction_analysis = reaction_analysis
        # set overall reaction analysis
        self.overall_reaction_analysis = overall_reaction_analysis

        # NOTE: set
        # component list
        self.component_list = list(component_dict.keys())

        # NOTE: kwargs
        self.threshold = kwargs.get('threshold', 1e-8)

    @property
    def eos(self):
        '''Get eos class'''
        # check
        return self._eos

    @eos.setter
    def eos(self, value: Any):
        '''Set eos class'''
        # set
        self._eos = value

    @property
    def eos_model(self):
        '''Get eos model'''
        return self._eos_model

    @eos_model.setter
    def eos_model(self, value: str):
        '''Set eos model'''
        # set
        self._eos_model = value

    @property
    def activity_model(self):
        '''Get activity model'''
        return self._activity_model

    @activity_model.setter
    def activity_model(self, value: str):
        '''Set activity model'''
        # set
        self._activity_model = value

    @property
    def activity(self):
        '''activity model (NRTL, UNIQUAC)'''
        return self._activity

    @activity.setter
    def activity(self, value: Any):
        '''Set activity model'''
        # set
        self._activity = value

    @property
    def activity_inputs(self):
        '''Get activity inputs'''
        return self._activity_inputs

    @activity_inputs.setter
    def activity_inputs(self, value: Dict[str, Any]):
        '''Set activity inputs'''
        # set
        self._activity_inputs = value

    @property
    def gas_mixture(self):
        '''Get gas mixture'''
        return self._gas_mixture

    @gas_mixture.setter
    def gas_mixture(self, value: str):
        '''Set gas mixture'''
        # set
        self._gas_mixture = value

    @property
    def solution(self):
        '''Get solution'''
        return self._solution

    @solution.setter
    def solution(self, value: str):
        '''Set solution'''
        # set
        self._solution = value

    def unpack_X(self, mol_data_pack: Dict[str, float | int]):
        '''
        convert X to dict {key: value} and list [value]

        Parameters
        ----------
        mol_data_pack : dict
            component mole dictionary {key: value}, such as {CO2: 0, H2: 1, CO: 2, H2O: 3, CH3OH: 4}

        Returns
        -------
        N0s_list : list
            The list of N0s
        N0s_vector : numpy.ndarray
            The vector of N0s
        Nf : float
            The total number of moles
        '''
        try:
            # size
            size = len(mol_data_pack)

            # lists
            N0s_list = []
            N0s_vector = np.zeros(size)

            # looping through
            for i, (com_key, com_value) in enumerate(self.component_dict.items()):
                N0s_list.append(mol_data_pack[com_key])
                N0s_vector[i] = mol_data_pack[com_key]

            # final
            Nf = np.sum(N0s_vector)

            return N0s_list, N0s_vector, Nf
        except Exception as e:
            raise Exception(
                f"Error in ReactionOptimizer.unpack_X(): {str(e)}") from e

    def build_EoR(self, EoR: List[float]):
        '''
        build EoR

        Parameters
        ----------
        EoR : list
            extent of reaction list

        Returns
        -------
        comp_value_list : list
            comp_value_list
        comp_value_matrix : numpy.ndarray
            comp_value_matrix
        '''
        try:
            # comp value list
            comp_value_list = []
            # new
            comp_list_updated = []

            # ! looping through reactions
            for i, item in enumerate(self.comp_list):
                _item = {}
                for key, value in item.items():
                    _item[key] = value * EoR[i]
                comp_list_updated.append(_item)

            # define a value list
            for i, item in enumerate(comp_list_updated):
                _value_list = []
                for key, value in item.items():
                    _value_list.append(value)
                comp_value_list.append(_value_list)

            # value matrix
            comp_value_matrix = np.array(comp_value_list)

            return comp_value_list, comp_value_matrix
        except Exception as e:
            raise Exception(
                f"Error in ReactionOptimizer.build_EoR(): {str(e)}") from e

    def build_final_X(
        self,
        N0s_vector: np.ndarray,
        EoR_vector: np.ndarray
    ):
        '''
        build final X

        Parameters
        ----------
        component_dict : dict
            component_dict
        N0s_vector : numpy.ndarray
            N0s_vector
        EoR_vector : numpy.ndarray
            EoR_vector

        Returns
        -------
        Xfs : dict
            Xfs
        Xfs_vector : numpy.ndarray
            Xfs_vector
        Nf : float
            Nf
        '''
        try:
            # final mole of components [mole]
            Nfs_vector = N0s_vector + EoR_vector

            # final mole fraction
            Xfs_vector = Nfs_vector/np.sum(Nfs_vector)

            # final mole
            Nf = np.sum(Nfs_vector)

            # convert to Xfs
            Xfs = {}
            for i, key in enumerate(self.component_dict.keys()):
                Xfs[key] = Xfs_vector[i]

            return Xfs, Xfs_vector, Nf, Nfs_vector
        except Exception as e:
            raise Exception(
                f"Error in ReactionOptimizer.build_final_X(): {str(e)}") from e

    def obj_fn(
        self,
        x,
        N0s: Dict[str, float],
        P: float,
        T: float,
        equilibrium_constants: Dict[str, float],
        method: Literal['minimize', 'least_squares'],
    ):
        '''
        Objective function for optimization

        Parameters
        ----------
        x : list
            loop values
        N0s : dict
            initial mole dictionary {key: value}, such as {CO2: 0, H2: 1, CO: 2, H2O: 3, CH3OH: 4}
            parameters
        P : float
            pressure [bar]
        T : float
            temperature [K]
        equilibrium_constants : dict
            reaction equilibrium constant dictionary calculated at T as: {key: value}, such as {R1: 0.1, R2: 0.2, ...}
        method : str
            optimization method.

        Returns
        -------
        obj : float | list
            objective function, if method is 'minimize' return float, if method is 'least_squares' return list
        '''
        try:
            # NOTE: default values
            # pressure [bar]
            # self.P_Ref_bar

            # NOTE: extent of reaction
            EoR = x

            # extent of reaction
            EoR = x
            # define threshold
            for i, item in enumerate(EoR):
                EoR[i] = max(self.threshold, item)

            # SECTION: mole balance
            # NOTE: initial mole
            # unpack N0s
            N0s_list, N0s_vector, N0f = self.unpack_X(N0s)

            # NOTE: update comp_list with EoR
            # build EoR => item[key] = value * EoR[i]
            _, comp_value_matrix = self.build_EoR(EoR)

            # NOTE: extent sun [mol]
            EoR_vector = np.sum(comp_value_matrix, axis=0)

            # NOTE: build final X
            Xfs, Xfs_vector, Nf, Nfs_vector = self.build_final_X(
                N0s_vector, EoR_vector)

            # SECTION: component state analysis

            # SECTION: equilibrium equation
            # reaction term
            reaction_terms = self.reaction_equilibrium_equation(
                Xfs=Xfs,
                P=P,
                T=T,
                equilibrium_constants=equilibrium_constants,
            )

            # SECTION: build objective functions
            if method == 'minimize':
                obj = 0
                for i, reaction in enumerate(self.reaction_analysis):
                    # NOTE: method 1
                    # obj = abs(reaction_terms[reaction])
                    # NOTE: method 2
                    # obj += reaction_terms[reaction]**2
                    # NOTE: set
                    obj += reaction_terms[reaction]**2

                # NOTE: define objective function
                # penalty obj
                # obj += 1e-3

                # sqrt
                # obj = sqrt(obj)
            elif method == 'least_squares':
                # NOTE: set
                obj = []
                for i, reaction in enumerate(self.reaction_analysis):
                    # NOTE: method 1
                    # obj.append(reaction_terms[reaction])
                    # NOTE: method 2
                    # obj.append(reaction_terms[reaction]**2)
                    # NOTE: set
                    obj.append(reaction_terms[reaction])

                # NOTE: mole fraction constraints
                Xfs_sum = np.sum(Xfs_vector)
                obj.append(Xfs_sum - 1)
            else:
                raise ValueError(
                    f"Invalid optimization method: {method}. Must be 'minimize' or 'least_squares'.")

            return obj
        except Exception as e:
            raise Exception(
                f"Error in objective function: {str(e)}") from e

    def reaction_equilibrium_equation(
        self,
        Xfs: Dict[str, float],
        P: float,
        T: float,
        equilibrium_constants: Dict[str, Dict[str, Any]],
        phase: Literal[
            "liquid", "gas"
        ] = "gas",
        **kwargs
    ):
        """
        Generate reaction equilibrium equations for gas and liquid phases.

        Parameters
        ----------
        Xfs : dict
            Mole fraction dictionary {key: value}, such as {CO2: 0, H2: 1, CO: 2, H2O: 3, CH3OH: 4}
        P : float
            Pressure [Pa]
        T : float
            Temperature [K]
        equilibrium_constants : dict
            Reaction equilibrium constant dictionary calculated at T as: such as {R1: dict, R2: dict, ...},
            the dict value consists of `value`, `symbol`, `unit`, `temperature`, `reaction`, `method`
        phase : str, optional
            Phase of calculation, by default "gas".
        kwargs : dict
            Additional parameters for non-ideal calculations.

        Notes
        -----
        The equilibrium constant at temperature T is calculated from the standard Gibbs free energy change:

        K_T = exp(-ΔG_T⁰ / RT) = ∏ [ (f̂_i / f_i⁰) ^ ν_i ]

        Where:
        - K_T    : Equilibrium constant at temperature T
        - ΔG_T⁰  : Standard Gibbs free energy change at temperature T (J/mol)
        - R      : Universal gas constant (8.314 J/mol·K)
        - T      : Temperature in Kelvin
        - f̂_i    : The mixture fugacity of component i.
        - f_i⁰   : The fugacity of pure component i in its reference state at temperature T.
        - ν_i    : Stoichiometric coefficient of component i (positive for products, negative for reactants)

        This expression connects thermodynamic data to fugacity-based equilibrium conditions.

        Gas Phase Calculation:
        - For gaseous compounds, the reference is the ideal gas at 1 bar, then f_i⁰ = P_ref = 1 bar.
        - f̂_i for non-ideal gases is calculated as f̂_i = φ_i * P * Y_i.
        - φ_i is the fugacity coefficient in the mixture.

        Liquid Phase Calculation:
        - f_i⁰ is standard-state fugacity (e.g., pure liquid at 1 bar)
        - f̂_i for liquid mixtures is defined as f̂_i = γ_i * X_i * f_i(pure,T).
        - f_i(pure,T) is the fugacity of pure component i at temperature T.
        - γ_i is the activity coefficient of component i in the liquid phase.
        - X_i is the mole fraction of component i in the liquid phase.

        Assumptions:
        - Regular liquid mixtures: Lewis-Randall
        - Ideal liquid solution: Raoult
        - Dilute or electrolyte systems: Henry
        - Poynting correction for liquid mixtures at high pressures.
        - General theoretical frameworks: Fugacity-based (chemical potential)
        """
        try:
            # NOTE: default values
            # pressure [bar]
            # self.P_Ref_bar
            # temperature [K]
            # self.T_Ref_K

            # SECTION: fugacity coefficient
            # fugacity coefficient
            fugacity_coeff = {}

            # check gas mixture
            if self.gas_mixture.lower() == "non-ideal":
                # NOTE: model input
                model_input = {
                    "feed-specification": Xfs,
                    "pressure": [P, 'bar'],
                    "temperature": [T, 'K'],
                }

                # NOTE: eos model
                eos_model = self.eos_model

                # Validate the model name
                if eos_model not in EOS_MODELS:
                    raise ValueError(
                        f"Invalid EOS model: {eos_model}. Must be {EOS_MODELS}.")

                # NOTE: calculate fugacity
                res_ = self._cal_fugacity_coefficient_gaseous_mixture(
                    model_name=eos_model,  # type: ignore
                    model_input=model_input)
                # update
                fugacity_coeff = {**res_}

            elif self.gas_mixture.lower() == "ideal":
                # set
                for i, key in enumerate(self.component_dict.keys()):
                    fugacity_coeff[key] = 1
            else:
                raise ValueError(
                    "Invalid gas mixture mode. Must be 'ideal' or 'non-ideal'.")

            # SECTION: activity coefficient
            # activity coefficient
            activity_coeff = {}

            # check liquid mixture
            if self.solution.lower() == "non-ideal":

                # prepare model input
                # Ensure activity_inputs is a dictionary before unpacking
                activity_inputs = self.activity_inputs or {}
                # set
                model_inputs = {
                    "mole_fraction": Xfs, **activity_inputs
                }

                # NOTE: check model name
                if self.activity_model == 'NRTL':
                    # exec
                    activity_coeff = self._cal_activity_coefficient_solution(
                        model_name='NRTL', model_input=model_inputs)
                elif self.activity_model == 'UNIQUAC':
                    # exec
                    activity_coeff = self._cal_activity_coefficient_solution(
                        model_name='UNIQUAC', model_input=model_inputs)
                else:
                    raise ValueError(
                        f"Invalid activity model: {self.activity_model}. Must be {ACTIVITY_MODELS}.")
            elif self.solution.lower() == "ideal":
                # set
                for i, key in enumerate(self.component_dict.keys()):
                    activity_coeff[key] = 1
            else:
                raise ValueError(
                    "Invalid liquid mixture mode. Must be 'ideal' or 'non-ideal'.")

            # SECTION: equilibrium equation
            # reaction term
            reaction_terms = {}

            # NOTE: loop over reactions
            for i, reaction in enumerate(self.reaction_analysis):
                # denominator
                denominator = 1
                # numerator
                numerator = 1

                # state count
                state_count_ = self.reaction_analysis[reaction]['state_count']

                # SECTION: loop over reactants
                for item in self.reaction_analysis[reaction]['reactants']:
                    # NOTE: item info
                    molecule_ = item['molecule']
                    molecule_state_ = item['molecule_state']
                    coefficient_ = item['coefficient']
                    state_ = item['state']

                    # final mole fraction
                    Xfs_ = Xfs[molecule_state_]

                    # NOTE: cal
                    # check phase
                    if state_ == "g":
                        # set
                        term_ = Xfs_ * (P / self.P_Ref_bar) * \
                            fugacity_coeff[molecule_state_]
                    elif state_ == "l":
                        # check solution
                        if self.solution.lower() == "non-ideal":
                            # Lewis-Randall/Raoult
                            term_ = Xfs_ * activity_coeff[molecule_state_]
                        elif self.solution.lower() == "ideal":
                            # check
                            if state_count_['l'] == 1:
                                # pure liquid
                                term_ = 1
                        else:
                            raise ValueError(
                                "Invalid liquid mixture mode. Must be 'ideal' or 'non-ideal'.")

                    elif state_ == "s":
                        # set
                        # solid always has the activity of 1
                        term_ = 1
                    else:
                        raise ValueError(
                            "Invalid phase. Must be 'gas' or 'liquid'.")

                    # update denominator
                    denominator *= (term_)**coefficient_

                # SECTION: loop over products
                for item in self.reaction_analysis[reaction]['products']:
                    # NOTE: item info
                    molecule_ = item['molecule']
                    molecule_state_ = item['molecule_state']
                    coefficient_ = item['coefficient']
                    state_ = item['state']

                    # final mole fraction
                    Xfs_ = Xfs[molecule_state_]

                    # NOTE: cal
                    # check phase
                    if state_ == "g":
                        # set
                        term_ = Xfs_ * (P / self.P_Ref_bar) * \
                            fugacity_coeff[molecule_state_]
                    elif state_ == "l":
                        # check solution
                        if self.solution.lower() == "non-ideal":
                            # Lewis-Randall/Raoult
                            term_ = Xfs_ * activity_coeff[molecule_state_]
                        elif self.solution.lower() == "ideal":
                            # check
                            if state_count_['l'] == 1:
                                # pure liquid
                                term_ = 1
                        else:
                            raise ValueError(
                                "Invalid liquid mixture mode. Must be 'ideal' or 'non-ideal'.")
                    elif state_ == "s":
                        # set
                        # solid always has the activity of 1
                        term_ = 1
                    else:
                        raise ValueError(
                            "Invalid phase. Must be 'gas' or 'liquid'.")

                    # update numerator
                    numerator *= (term_)**coefficient_

                # NOTE: update reaction term
                # equilibrium data
                Keq_ = equilibrium_constants[reaction]['value']

                # ? method 1
                # reaction_terms[reaction] = (numerator/denominator) - Keq_
                # ? method 2
                # reaction_terms[reaction] = numerator - denominator*Keq_
                # ? method 3
                Q_i = numerator / denominator
                Q_i_safe = np.clip(Q_i, 1e-12, None)
                reaction_terms[reaction] = log(Q_i_safe) - log(Keq_)

            # return
            return reaction_terms
        except Exception as e:
            raise Exception(
                f"Error in building the reaction equilibrium equation for {phase}: {str(e)}") from e

    def _cal_fugacity_coefficient_gaseous_mixture(
        self,
        model_name: Literal[
            "SRK", "PR", "RK"
        ],
        model_input: Dict[str, Any]
    ):
        """
        Calculate the fugacity coefficient of gaseous mixture using the specified EOS model.

        Parameters
        ----------
        model_name : str
            The EOS model to use for calculation. Options: "SRK", "PR", "RK".
        model_input : dict
            The input data for the EOS model.
            - feed-specification: Dictionary of component mole fractions.
            - pressure: pressure value in any unit.
            - temperature: temperature value in any unit.

        Returns
        -------
        dict
            A dictionary containing the fugacity coefficients for each component in the gaseous mixture.
        """
        try:
            # SECTION: model source
            model_source = {
                "datasource": self.datasource,
                "equationsource": self.equationsource
            }

            # NOTE: calculate fugacity
            eos = self.eos

            # check
            if eos is None or eos == 'None':
                raise ValueError(
                    f"Invalid EOS model. Must be {EOS_MODELS}.")

            res = eos.cal_fugacity_mixture(
                model_name=model_name,
                model_input=model_input,
                model_source=model_source
            )

            # NOTE: extract fugacity coefficient for each component
            res_1 = res['vapor']

            # res
            phi_comp = {}

            # looping through each component
            for i, key in enumerate(res_1.keys()):
                # set
                phi_comp[key] = res_1[key]['fugacity_coefficient']['value']

            return phi_comp
        except Exception as e:
            raise Exception(
                f"Error in calculating the fugacity coefficient for the gaseous mixture: {str(e)}") from e

    def _cal_activity_coefficient_solution(
        self,
        model_name: Literal[
            'NRTL', 'UNIQUAC'
        ],
        model_input: Dict[str, Any],
    ):
        """
        Calculate the activity coefficient of solution using the specified model.

        Parameters
        ----------
        model_name : str
            The model to use for calculation. Options: "NRTL", "UNIQUAC".
        model_input : dict
            The input data for the model.
            - mole_fraction: Dictionary of component mole fractions.
            - tau_ij: Dictionary of interaction parameters.
            - alpha_ij: Dictionary of interaction parameters.
            - r_i: relative van der Waals volume of component i
            - q_i: relative surface area of component i
        model_source : dict
            The source of the model data.
            - datasource: Data source for the model.
            - equationsource: Equation source for the model.

        Returns
        -------
        dict
            A dictionary containing the fugacity coefficients for each component in the liquid mixture.
        """
        try:
            # SECTION: model source
            # model_source = {
            #     "datasource": self.datasource,
            #     "equationsource": self.equationsource
            # }
            # check activity
            if self.activity is None or self.activity == 'None':
                raise ValueError(
                    f"Invalid activity model. Must be {ACTIVITY_MODELS}.")

            # NOTE: model name
            if model_name == "NRTL":
                # calculate fugacity
                res, others = self.activity.cal(model_input=model_input)
            elif model_name == "UNIQUAC":
                # calculate fugacity
                res, others = self.activity.cal(model_input=model_input)
            else:
                raise ValueError(
                    f"Invalid model name: {model_name}. Must be 'NRTL' or 'UNIQUAC'.")

            # NOTE: extract fugacity coefficient for each component
            return others['AcCo_i_comp']
        except Exception as e:
            raise Exception(
                f"Error in calculating the fugacity coefficient for the liquid mixture: {str(e)}") from e

    def constraint1(self, x, params):
        '''
        Sets a bound between 0 and 1 [x>>0 & x<<1] representing {-f(x) + 1 >= 0}

        Parameters
        ------------
        x : list
            loop values
        params : list
            parameters

        Returns
        -------
        obj : float
            objective function

        Notes
        -----
        1. component mole fraction should be less than 1
        '''
        # print(f"1: {x}")
        N0s, comp_list, component_dict, i = params

        # build EoR
        comp_value_list, comp_value_matrix = self.build_EoR(x)

        # extent sun [mol]
        EoR_vector = np.sum(comp_value_matrix, axis=0)

        # unpack N0s
        N0s_list, N0s_vector, N0f = self.unpack_X(N0s)

        # build final X
        Xfs, Xfs_vector, Nf, Nfs_vector = self.build_final_X(
            N0s_vector, EoR_vector)

        # constraint
        cons = -1*Xfs_vector[i] + 1

        return cons

    def constraint2(self, x, params):
        '''
        Sets a bound between 0 and 1 [x>>0 & x<<1] representing {f(x) >= 0}

        Parameters
        ------------
        x : list
        loop values
        params : list
        parameters

        Returns
        -------
        obj : float
        objective function

        Notes
        -----
        1. component mole fraction should be greater than 0
        2. define epsilon constraint 1e-5
        '''
        # print(f"2: {x}")
        N0s, comp_list, component_dict, i = params

        # build EoR
        comp_value_list, comp_value_matrix = self.build_EoR(x)

        # extent sun [mol]
        EoR_vector = np.sum(comp_value_matrix, axis=0)

        # unpack N0s
        N0s_list, N0s_vector, N0f = self.unpack_X(N0s)

        # build final X
        Xfs, Xfs_vector, Nf, Nfs_vector = self.build_final_X(
            N0s_vector, EoR_vector)

        # constraint
        # ! set epsilon constraint
        cos = Xfs_vector[i] - 1e-5

        return cos

    def constraint3(self, x, params):
        '''
        Sets a bound for total mole f(x)>=0

        Parameters
        ------------
        x : list
        loop values
        params : list
        parameters

        Returns
        -------
        cons : float
        constraint

        Notes
        -----
        1. total mole of components should be greater than 0
        '''
        # print(f"3: {x}")
        N0s, comp_list, component_dict, i = params

        # build EoR
        _, comp_value_matrix = self.build_EoR(x)

        # extent sun [mol]
        EoR_vector = np.sum(comp_value_matrix, axis=0)

        # unpack N0s
        N0s_list, N0s_vector, N0f = self.unpack_X(N0s)

        # build final X
        Xfs, Xfs_vector, Nf, Nfs_vector = self.build_final_X(
            N0s_vector, EoR_vector)

        # cons
        cons = Nfs_vector[i]

        return cons

    def constraint4(self, x, params):
        '''
        Set a bound for reactants which are consumped ff(x)<=f0(x)

        Parameters
        ------------
        x : list
        loop values
        params : list
        parameters

        Returns
        -------
        cons : float
        constraint

        Notes
        -----
        1. final mole of reactants should be greater than 0
        '''
        # print(f"4: {x}")
        N0s, comp_list, component_dict, i = params

        # build EoR
        comp_value_list, comp_value_matrix = self.build_EoR(x)

        # extent sun [mol]
        EoR_vector = np.sum(comp_value_matrix, axis=0)

        # unpack N0s
        N0s_list, N0s_vector, N0f = self.unpack_X(N0s)

        # build final X
        Xfs, Xfs_vector, Nf, Nfs_vector = self.build_final_X(
            N0s_vector, EoR_vector)

        # cons
        cons = -1*Nfs_vector[i] + N0s_vector[i]

        return cons

    def constraint5(self, x, params):
        '''
        Sets a bound for sum of all mole fraction (x[i]<1)

        Parameters
        ------------
        x : list
        loop values
        params : list
        parameters

        Returns
        -------
        cons : float
        constraint

        Notes
        -----
        1. sum of all mole fraction should be 1
        '''
        # print(f"5: {x}")
        N0s, comp_list, component_dict = params

        # build EoR
        comp_value_list, comp_value_matrix = self.build_EoR(x)

        # extent sun [mol]
        EoR_vector = np.sum(comp_value_matrix, axis=0)

        # unpack N0s
        N0s_list, N0s_vector, N0f = self.unpack_X(N0s)

        # build final X
        Xfs, Xfs_vector, Nf, Nfs_vector = self.build_final_X(
            N0s_vector, EoR_vector)

        # cons
        cons = np.sum(Xfs_vector) - 1

        return cons

    def constraint6(self, x, params):
        '''
        Sets a bound for sum of all mole fraction (x[i] = 1) representing f(x)=1

        Parameters
        ------------
        x : list
        loop values
        params : list
        parameters

        Returns
        -------
        cons : float
        constraint

        Notes
        -----
        1. sum of all mole fraction should be 1
        '''
        # print(f"5: {x}")
        N0s, comp_list, component_dict = params

        # build EoR
        comp_value_list, comp_value_matrix = self.build_EoR(x)

        # extent sun [mol]
        EoR_vector = np.sum(comp_value_matrix, axis=0)

        # unpack N0s
        N0s_list, N0s_vector, N0f = self.unpack_X(N0s)

        # build final X
        Xfs, Xfs_vector, Nf, Nfs_vector = self.build_final_X(
            N0s_vector, EoR_vector)

        # cons
        cons = np.sum(Xfs_vector) - 1

        return cons

    def opt_run(
        self,
        initial_mole: Dict[str, float | int],
        initial_mole_fraction: Dict[str, float | int],
        pressure: float,
        temperature: float,
        equilibrium_constants: Dict[str, float],
        reaction_numbers: int,
        method: Literal['minimize', 'least_squares'] = 'minimize',
        **kwargs
    ):
        """
        Start optimization process

        Parameters
        ----------
        initial_mole : dict
            Initial mole dictionary {key: value}, such as {CO2: 0, H2: 1, CO: 2, H2O: 3, CH3OH: 4}
        initial_mole_fraction : dict
            Initial mole fraction dictionary {key: value}, such as {CO2: 0, H2: 1, CO: 2, H2O: 3, CH3OH: 4}
        pressure : float
            Pressure [bar]
        temperature : float
            Temperature [K]
        equilibrium_constants : dict
            Reaction equilibrium constant dictionary calculated at T as: {key: value}, such as {R1: 0.1, R2: 0.2, ...}
        reaction_numbers : int
            Number of reactions
        method : str, optional
            Optimization method, by default 'minimize'.
        kwargs : dict
            Additional parameters for optimization.
            - least_square_algorithm: 'trf', 'dogbox', 'lm'
            - minimize_algorithm: 'SLSQP', 'trust-constr', 'COBYLA'
            - scale: Scale factor for the upper bound.
            - bound_scale: Scale factor for the upper bound.
            - fallback: Fallback value for upper bound if no reactants are present.

        Returns
        -------
        opt_res : OptimizeResult
            Optimization result object containing the optimized values and other information.
        """
        try:
            # NOTE: default values
            # least_square_algorithm
            least_square_algorithm = kwargs.get(
                'least_square_algorithm', 'trf')
            # minimize_algorithm
            minimize_algorithm = kwargs.get('minimize_algorithm', 'SLSQP')

            # ? check if the algorithm is valid
            if minimize_algorithm not in ['SLSQP', 'trust-constr', 'COBYLA']:
                raise ValueError(
                    f"Invalid minimize algorithm: {minimize_algorithm}. Must be 'SLSQP', 'trust-constr', or 'COBYLA'.")

            # ? check if the algorithm is valid
            if least_square_algorithm not in ['trf', 'dogbox', 'lm']:
                raise ValueError(
                    f"Invalid least square algorithm: {least_square_algorithm}. Must be 'trf', 'dogbox', or 'lm'.")

            # NOTE: set
            # initial guess for extent of reaction
            # default value
            EOR0_val = kwargs.get('EOR0_val', 0.5)

            # set
            EOR0 = np.random.uniform(0, EOR0_val, reaction_numbers)

            # NOTE: bounds
            bound0 = (0, 20)
            bounds = []
            for i in range(len(EOR0)):
                bounds.append(bound0)

            # EOR02
            # initial mole
            initial_mole_ = np.array(list(initial_mole.values()))

            # extent of reaction (EoR)
            EOR0_Bounds, EOR0_bounds, EoR_initial = self.compute_bounds(
                nu=self.stoichiometric_coeff,
                n0=initial_mole_,
                **kwargs,
            )

            # NOTE: constraints
            # constraints 1
            cons1 = self.constraints_collection_1(initial_mole)
            # constraints 2
            cons2 = self.constraints_collection_2(
                nu=self.stoichiometric_coeff,
                n0=initial_mole_,
                include_mole_fraction_constraint=True
            )

            # set constraints
            cons = []

            # SECTION: optimize
            if method == 'minimize':
                opt_res = optimize.minimize(
                    fun=self.obj_fn,
                    x0=EOR0,
                    args=(
                        initial_mole,
                        pressure,
                        temperature,
                        equilibrium_constants,
                        method
                    ),
                    method=minimize_algorithm,
                    bounds=EOR0_bounds,
                    constraints=cons1,
                    options={
                        'disp': False,
                        'ftol': 1e-12,
                        'maxiter': 1000
                    }
                )
            elif method == 'least_squares':
                opt_res = optimize.least_squares(
                    fun=self.obj_fn,
                    x0=EOR0,
                    args=(
                        initial_mole,
                        pressure,
                        temperature,
                        equilibrium_constants,
                        method
                    ),
                    bounds=EOR0_Bounds,
                    method=least_square_algorithm,
                )
            else:
                raise ValueError(
                    f"Invalid optimization method: {method}. Must be 'minimize' or 'least_squares'.")

            # save
            return opt_res
        except Exception as e:
            raise Exception(
                f"Error in the optimization process: {str(e)}") from e

    def constraints_collection_1(self, initial_mole: Dict[str, float | int]):
        """
        Generate a collection of constraints for optimization.

        Parameters
        ----------
        initial_mole : Dict[str, float | int]
            Initial mole dictionary {key: value}, such as {CO2: 0, H2: 1, CO: 2, H2O: 3, CH3OH: 4}
        """
        try:
            # NOTE: define constraint
            cons = []

            # inequality
            # looping through components
            for key, value in self.component_dict.items():
                # append constraint
                cons.append({
                    'type': 'ineq',
                    'fun': self.constraint1,
                    'args': (
                        (initial_mole, self.comp_list, self.component_dict, value),)
                })
                # append constraint
                cons.append({
                    'type': 'ineq',
                    'fun': self.constraint2,
                    'args': (
                        (initial_mole, self.comp_list, self.component_dict, value),)
                })
                # append constraint
                cons.append({
                    'type': 'ineq',
                    'fun': self.constraint3,
                    'args': (
                        (initial_mole, self.comp_list, self.component_dict, value),)
                })

                # check consumed
                if key in self.overall_reaction_analysis['consumed']:
                    # append constraint
                    cons.append({
                        'type': 'ineq',
                        'fun': self.constraint4,
                        'args': (
                            (initial_mole, self.comp_list, self.component_dict, value),)
                    })

            # append constraint
            cons.append({
                'type': 'ineq',
                'fun': self.constraint5,
                        'args': ((initial_mole, self.comp_list, self.component_dict),)
            })

            # equality
            # constraint 5
            # cons.append({
            # 'type': 'ineq',
            # 'fun': constraint5,
            # 'args':((N0s, comp_list,component_dict),)
            # })

            # return
            return cons
        except Exception as e:
            raise Exception(
                f"Error in generating constraints collection: {str(e)}") from e

    def constraints_collection_2(
        self,
        nu,
        n0,
        include_mole_fraction_constraint=True
    ):
        """
        Generate a collection of constraints for optimization.

        Parameters
        ----------
        initial_mole : Dict[str, float | int]
            Initial mole dictionary {key: value}, such as {CO2: 0, H2: 1, CO: 2, H2O: 3, CH3OH: 4}
        """
        try:
            n_species, n_reactions = nu.shape
            constraints = []

            # 1. Species non-negativity constraints
            for i in range(n_species):
                def make_fn(i):
                    return lambda ksi: n0[i] + np.dot(nu[i, :], ksi)
                constraints.append(NonlinearConstraint(make_fn(i), 0, np.inf))

            # 2. Mole fraction sum constraint (optional)
            if include_mole_fraction_constraint:
                def mole_fraction_sum(ksi):
                    n = n0 + nu @ ksi
                    total = np.sum(n)
                    # avoid division by zero if total is very small
                    if total <= 1e-12:
                        return 0.0
                    y = n / total
                    return np.sum(y)

                constraints.append(NonlinearConstraint(
                    mole_fraction_sum, 1.0, 1.0))

            return constraints
        except Exception as e:
            raise Exception(
                f"Error in generating constraints collection: {str(e)}") from e

    def equilibrium_results(
        self,
        initial_mole: Dict[str, float | int],
        EoR: list
    ):
        '''
        Calculate the equilibrium results based on the initial mole and extent of reaction.

        Parameters
        ----------
        initial_mole : dict
            Initial mole dictionary {key: value}, such as {CO2: 0, H2: 1, CO: 2, H2O: 3, CH3OH: 4}
        EoR : list
            Extent of reaction list [EoR1, EoR2, EoR3, ...]

        Returns
        -------
        Xfs : dict
            Final mole fraction dictionary {key: value}, such as {CO2: 0.1, H2: 0.1, CO: 0.2, H2O: 0.3, CH3OH: 0.4}
        Xfs_vector : np.ndarray
            Final mole fraction vector [Xfs1, Xfs2, Xfs3, ...]
        Nf : float
            Final moles of the system
        Nfs_vector : np.ndarray
            Final moles vector [Nf1, Nf2, Nf3, ...]
        Nfs : dict
            Final moles dictionary {key: value}, such as {CO2: 0, H2: 1, CO: 2, H2O: 3, CH3OH: 4}
        '''
        try:
            # unpack
            N0s_list, N0s_vector, N0f = self.unpack_X(initial_mole)

            # update comp_list with EoR
            # build EoR => item[key] = value * EoR[i]
            _, comp_value_matrix = self.build_EoR(EoR)

            # extent sun [mol]
            EoR_vector = np.sum(comp_value_matrix, axis=0)

            # build final X
            Xfs, Xfs_vector, Nf, Nfs_vector = self.build_final_X(
                N0s_vector, EoR_vector
            )

            # Nfs
            Nfs = {}
            for key, value in self.component_dict.items():
                Nfs[key] = Nfs_vector[value]

            # set float
            Nfs = {k: float(v) for k, v in Nfs.items()}
            Xfs = {k: float(v) for k, v in Xfs.items()}

            return Xfs, Xfs_vector, Nf, Nfs_vector, Nfs
        except Exception as e:
            raise Exception(
                f"Error in processing optimization results: {str(e)}") from e

    def compute_bounds(
        self,
        nu: np.ndarray,
        n0: np.ndarray,
        **kwargs
    ):
        """
        Computes lower and upper bounds for each reaction extent ξ_j.

        Parameters
        ----------
        nu: np.ndarray
            The stoichiometric matrix (n_species x n_reactions)
        N0s: np.ndarray
            Initial moles of each species (length: n_species)
        kwargs: dict
            Additional parameters for bounds calculation.
            - scale: float
                Scale factor for the upper bound.
            - bound_scale: float
                Scale factor for the upper bound.
            - fallback: float
                Fallback value for upper bound if no reactants are present.

        Returns
        -------
        Bounds
            scipy.optimize.Bounds object for use in minimize/least_squares.
        """
        try:
            # SECTION: default values
            # scale
            scale = kwargs.get('scale', 10.0)
            # bound_scale
            bound_scale = kwargs.get('bound_scale', 0.5)
            # fallback
            fallback = kwargs.get('fallback', None)

            # SECTION: Check if nu is a 2D array
            if nu.ndim != 2:
                raise ValueError("nu must be a 2D array.")
            # Check if n0 is a 1D array
            if n0.ndim != 1:
                raise ValueError("n0 must be a 1D array.")

            # NOTE: Check if nu and n0 have compatible dimensions
            if nu.shape[0] != n0.shape[0]:
                raise ValueError("nu and n0 must have compatible dimensions.")

            # NOTE: set fallback
            if fallback is None:
                fallback = max(n0.sum(), 1.0) * scale  # 10× total initial mol

            # SECTION: Compute bounds
            n_species, n_reactions = nu.shape
            lb = np.zeros(n_reactions)
            ub = np.full(n_reactions, np.inf)

            # Loop over reactions
            for j in range(n_reactions):
                for i in range(n_species):
                    v = nu[i, j]
                    if v < 0:
                        if n0[i] > 0:
                            max_mu = n0[i] / abs(v)
                            ub[j] = min(ub[j], max_mu)
                        else:
                            # No initial amount of a required reactant → can't proceed
                            ub[j] = min(ub[j], fallback)

            # NOTE: convert to bounds list
            bounds = []
            for i in range(len(lb)):
                # add slack to the upper bound
                ub[i] = ub[i]
                bounds.append((lb[i], ub[i]))

            # NOTE: initial guess
            EOR0 = []
            for i in range(len(ub)):
                # check if upper bound is finite
                if ub[i] == scale:
                    # set
                    _val = lb[i] + bound_scale * (ub[i] - lb[i])
                    EOR0.append(_val)
                else:
                    # set
                    EOR0.append(np.random.uniform(lb[i], ub[i]))

            return Bounds(lb, ub), bounds, EOR0
        except Exception as e:
            raise Exception(
                f"Error in computing bounds: {str(e)}") from e

activity property writable

activity model (NRTL, UNIQUAC)

activity_inputs property writable

Get activity inputs

activity_model property writable

Get activity model

eos property writable

Get eos class

eos_model property writable

Get eos model

gas_mixture property writable

Get gas mixture

solution property writable

Get solution

__init__(datasource, equationsource, component_dict, comp_list, stoichiometric_coeff, reaction_analysis, overall_reaction_analysis, **kwargs)

Initialize ReactionOptimizer class

Parameters

datasource : dict datasource dictionary equationsource : dict equationsource dictionary component_dict : dict component dictionary {key: value}, such as {CO2: 0, H2: 1, CO: 2, H2O: 3, CH3OH: 4} comp_list : list component list [{key: value}, {key: value}, ...], such as [{CO2: -1.0, H2: -3.0, CH3OH: 1.0, H2O: 1.0}, ...] reaction_analysis : dict reaction analysis result overall_reaction_analysis : dict overall reaction analysis result kwargs : dict additional parameters

Source code in PyReactLab/docs/optim.py

def __init__(
    self,
    datasource: Dict[str, Any],
    equationsource: Dict[str, Any],
    component_dict: Dict[str, float | int],
    comp_list: List[Dict[str, float | int]],
    stoichiometric_coeff: np.ndarray,
    reaction_analysis: Dict,
    overall_reaction_analysis: Dict,
    **kwargs
):
    '''
    Initialize ReactionOptimizer class

    Parameters
    ----------
    datasource : dict
        datasource dictionary
    equationsource : dict
        equationsource dictionary
    component_dict : dict
        component dictionary {key: value}, such as {CO2: 0, H2: 1, CO: 2, H2O: 3, CH3OH: 4}
    comp_list : list
        component list [{key: value}, {key: value}, ...], such as [{CO2: -1.0, H2: -3.0, CH3OH: 1.0, H2O: 1.0}, ...]
    reaction_analysis : dict
        reaction analysis result
    overall_reaction_analysis : dict
        overall reaction analysis result
    kwargs : dict
        additional parameters
    '''
    # datasource
    self.datasource = datasource
    # equationsource
    self.equationsource = equationsource
    # set component dictionary
    self.component_dict = component_dict
    # set component list
    self.comp_list = comp_list
    # set stoichiometric coefficient
    self.stoichiometric_coeff = stoichiometric_coeff
    # set reaction analysis
    self.reaction_analysis = reaction_analysis
    # set overall reaction analysis
    self.overall_reaction_analysis = overall_reaction_analysis

    # NOTE: set
    # component list
    self.component_list = list(component_dict.keys())

    # NOTE: kwargs
    self.threshold = kwargs.get('threshold', 1e-8)

build_EoR(EoR)

build EoR

Parameters

EoR : list extent of reaction list

Returns

comp_value_list : list comp_value_list comp_value_matrix : numpy.ndarray comp_value_matrix

Source code in PyReactLab/docs/optim.py

def build_EoR(self, EoR: List[float]):
    '''
    build EoR

    Parameters
    ----------
    EoR : list
        extent of reaction list

    Returns
    -------
    comp_value_list : list
        comp_value_list
    comp_value_matrix : numpy.ndarray
        comp_value_matrix
    '''
    try:
        # comp value list
        comp_value_list = []
        # new
        comp_list_updated = []

        # ! looping through reactions
        for i, item in enumerate(self.comp_list):
            _item = {}
            for key, value in item.items():
                _item[key] = value * EoR[i]
            comp_list_updated.append(_item)

        # define a value list
        for i, item in enumerate(comp_list_updated):
            _value_list = []
            for key, value in item.items():
                _value_list.append(value)
            comp_value_list.append(_value_list)

        # value matrix
        comp_value_matrix = np.array(comp_value_list)

        return comp_value_list, comp_value_matrix
    except Exception as e:
        raise Exception(
            f"Error in ReactionOptimizer.build_EoR(): {str(e)}") from e

build_final_X(N0s_vector, EoR_vector)

build final X

Parameters

component_dict : dict component_dict N0s_vector : numpy.ndarray N0s_vector EoR_vector : numpy.ndarray EoR_vector

Returns

Xfs : dict Xfs Xfs_vector : numpy.ndarray Xfs_vector Nf : float Nf

Source code in PyReactLab/docs/optim.py

def build_final_X(
    self,
    N0s_vector: np.ndarray,
    EoR_vector: np.ndarray
):
    '''
    build final X

    Parameters
    ----------
    component_dict : dict
        component_dict
    N0s_vector : numpy.ndarray
        N0s_vector
    EoR_vector : numpy.ndarray
        EoR_vector

    Returns
    -------
    Xfs : dict
        Xfs
    Xfs_vector : numpy.ndarray
        Xfs_vector
    Nf : float
        Nf
    '''
    try:
        # final mole of components [mole]
        Nfs_vector = N0s_vector + EoR_vector

        # final mole fraction
        Xfs_vector = Nfs_vector/np.sum(Nfs_vector)

        # final mole
        Nf = np.sum(Nfs_vector)

        # convert to Xfs
        Xfs = {}
        for i, key in enumerate(self.component_dict.keys()):
            Xfs[key] = Xfs_vector[i]

        return Xfs, Xfs_vector, Nf, Nfs_vector
    except Exception as e:
        raise Exception(
            f"Error in ReactionOptimizer.build_final_X(): {str(e)}") from e

compute_bounds(nu, n0, **kwargs)

Computes lower and upper bounds for each reaction extent ξ_j.

Parameters

nu: np.ndarray The stoichiometric matrix (n_species x n_reactions) N0s: np.ndarray Initial moles of each species (length: n_species) kwargs: dict Additional parameters for bounds calculation. - scale: float Scale factor for the upper bound. - bound_scale: float Scale factor for the upper bound. - fallback: float Fallback value for upper bound if no reactants are present.

Returns

Bounds scipy.optimize.Bounds object for use in minimize/least_squares.

Source code in PyReactLab/docs/optim.py

def compute_bounds(
    self,
    nu: np.ndarray,
    n0: np.ndarray,
    **kwargs
):
    """
    Computes lower and upper bounds for each reaction extent ξ_j.

    Parameters
    ----------
    nu: np.ndarray
        The stoichiometric matrix (n_species x n_reactions)
    N0s: np.ndarray
        Initial moles of each species (length: n_species)
    kwargs: dict
        Additional parameters for bounds calculation.
        - scale: float
            Scale factor for the upper bound.
        - bound_scale: float
            Scale factor for the upper bound.
        - fallback: float
            Fallback value for upper bound if no reactants are present.

    Returns
    -------
    Bounds
        scipy.optimize.Bounds object for use in minimize/least_squares.
    """
    try:
        # SECTION: default values
        # scale
        scale = kwargs.get('scale', 10.0)
        # bound_scale
        bound_scale = kwargs.get('bound_scale', 0.5)
        # fallback
        fallback = kwargs.get('fallback', None)

        # SECTION: Check if nu is a 2D array
        if nu.ndim != 2:
            raise ValueError("nu must be a 2D array.")
        # Check if n0 is a 1D array
        if n0.ndim != 1:
            raise ValueError("n0 must be a 1D array.")

        # NOTE: Check if nu and n0 have compatible dimensions
        if nu.shape[0] != n0.shape[0]:
            raise ValueError("nu and n0 must have compatible dimensions.")

        # NOTE: set fallback
        if fallback is None:
            fallback = max(n0.sum(), 1.0) * scale  # 10× total initial mol

        # SECTION: Compute bounds
        n_species, n_reactions = nu.shape
        lb = np.zeros(n_reactions)
        ub = np.full(n_reactions, np.inf)

        # Loop over reactions
        for j in range(n_reactions):
            for i in range(n_species):
                v = nu[i, j]
                if v < 0:
                    if n0[i] > 0:
                        max_mu = n0[i] / abs(v)
                        ub[j] = min(ub[j], max_mu)
                    else:
                        # No initial amount of a required reactant → can't proceed
                        ub[j] = min(ub[j], fallback)

        # NOTE: convert to bounds list
        bounds = []
        for i in range(len(lb)):
            # add slack to the upper bound
            ub[i] = ub[i]
            bounds.append((lb[i], ub[i]))

        # NOTE: initial guess
        EOR0 = []
        for i in range(len(ub)):
            # check if upper bound is finite
            if ub[i] == scale:
                # set
                _val = lb[i] + bound_scale * (ub[i] - lb[i])
                EOR0.append(_val)
            else:
                # set
                EOR0.append(np.random.uniform(lb[i], ub[i]))

        return Bounds(lb, ub), bounds, EOR0
    except Exception as e:
        raise Exception(
            f"Error in computing bounds: {str(e)}") from e

constraint1(x, params)

Sets a bound between 0 and 1 [x>>0 & x<<1] representing {-f(x) + 1 >= 0}

Parameters

x : list loop values params : list parameters

Returns

obj : float objective function

Notes

1. component mole fraction should be less than 1

Source code in PyReactLab/docs/optim.py

def constraint1(self, x, params):
    '''
    Sets a bound between 0 and 1 [x>>0 & x<<1] representing {-f(x) + 1 >= 0}

    Parameters
    ------------
    x : list
        loop values
    params : list
        parameters

    Returns
    -------
    obj : float
        objective function

    Notes
    -----
    1. component mole fraction should be less than 1
    '''
    # print(f"1: {x}")
    N0s, comp_list, component_dict, i = params

    # build EoR
    comp_value_list, comp_value_matrix = self.build_EoR(x)

    # extent sun [mol]
    EoR_vector = np.sum(comp_value_matrix, axis=0)

    # unpack N0s
    N0s_list, N0s_vector, N0f = self.unpack_X(N0s)

    # build final X
    Xfs, Xfs_vector, Nf, Nfs_vector = self.build_final_X(
        N0s_vector, EoR_vector)

    # constraint
    cons = -1*Xfs_vector[i] + 1

    return cons

constraint2(x, params)

Sets a bound between 0 and 1 [x>>0 & x<<1] representing {f(x) >= 0}

Parameters

x : list loop values params : list parameters

Returns

obj : float objective function

Notes

1. component mole fraction should be greater than 0
2. define epsilon constraint 1e-5

Source code in PyReactLab/docs/optim.py

def constraint2(self, x, params):
    '''
    Sets a bound between 0 and 1 [x>>0 & x<<1] representing {f(x) >= 0}

    Parameters
    ------------
    x : list
    loop values
    params : list
    parameters

    Returns
    -------
    obj : float
    objective function

    Notes
    -----
    1. component mole fraction should be greater than 0
    2. define epsilon constraint 1e-5
    '''
    # print(f"2: {x}")
    N0s, comp_list, component_dict, i = params

    # build EoR
    comp_value_list, comp_value_matrix = self.build_EoR(x)

    # extent sun [mol]
    EoR_vector = np.sum(comp_value_matrix, axis=0)

    # unpack N0s
    N0s_list, N0s_vector, N0f = self.unpack_X(N0s)

    # build final X
    Xfs, Xfs_vector, Nf, Nfs_vector = self.build_final_X(
        N0s_vector, EoR_vector)

    # constraint
    # ! set epsilon constraint
    cos = Xfs_vector[i] - 1e-5

    return cos

constraint3(x, params)

Sets a bound for total mole f(x)>=0

Parameters

x : list loop values params : list parameters

Returns

cons : float constraint

Notes

1. total mole of components should be greater than 0

Source code in PyReactLab/docs/optim.py

def constraint3(self, x, params):
    '''
    Sets a bound for total mole f(x)>=0

    Parameters
    ------------
    x : list
    loop values
    params : list
    parameters

    Returns
    -------
    cons : float
    constraint

    Notes
    -----
    1. total mole of components should be greater than 0
    '''
    # print(f"3: {x}")
    N0s, comp_list, component_dict, i = params

    # build EoR
    _, comp_value_matrix = self.build_EoR(x)

    # extent sun [mol]
    EoR_vector = np.sum(comp_value_matrix, axis=0)

    # unpack N0s
    N0s_list, N0s_vector, N0f = self.unpack_X(N0s)

    # build final X
    Xfs, Xfs_vector, Nf, Nfs_vector = self.build_final_X(
        N0s_vector, EoR_vector)

    # cons
    cons = Nfs_vector[i]

    return cons

constraint4(x, params)

Set a bound for reactants which are consumped ff(x)<=f0(x)

Parameters

x : list loop values params : list parameters

Returns

cons : float constraint

Notes

1. final mole of reactants should be greater than 0

Source code in PyReactLab/docs/optim.py

def constraint4(self, x, params):
    '''
    Set a bound for reactants which are consumped ff(x)<=f0(x)

    Parameters
    ------------
    x : list
    loop values
    params : list
    parameters

    Returns
    -------
    cons : float
    constraint

    Notes
    -----
    1. final mole of reactants should be greater than 0
    '''
    # print(f"4: {x}")
    N0s, comp_list, component_dict, i = params

    # build EoR
    comp_value_list, comp_value_matrix = self.build_EoR(x)

    # extent sun [mol]
    EoR_vector = np.sum(comp_value_matrix, axis=0)

    # unpack N0s
    N0s_list, N0s_vector, N0f = self.unpack_X(N0s)

    # build final X
    Xfs, Xfs_vector, Nf, Nfs_vector = self.build_final_X(
        N0s_vector, EoR_vector)

    # cons
    cons = -1*Nfs_vector[i] + N0s_vector[i]

    return cons

constraint5(x, params)

Sets a bound for sum of all mole fraction (x[i]<1)

Parameters

x : list loop values params : list parameters

Returns

cons : float constraint

Notes

1. sum of all mole fraction should be 1

Source code in PyReactLab/docs/optim.py

def constraint5(self, x, params):
    '''
    Sets a bound for sum of all mole fraction (x[i]<1)

    Parameters
    ------------
    x : list
    loop values
    params : list
    parameters

    Returns
    -------
    cons : float
    constraint

    Notes
    -----
    1. sum of all mole fraction should be 1
    '''
    # print(f"5: {x}")
    N0s, comp_list, component_dict = params

    # build EoR
    comp_value_list, comp_value_matrix = self.build_EoR(x)

    # extent sun [mol]
    EoR_vector = np.sum(comp_value_matrix, axis=0)

    # unpack N0s
    N0s_list, N0s_vector, N0f = self.unpack_X(N0s)

    # build final X
    Xfs, Xfs_vector, Nf, Nfs_vector = self.build_final_X(
        N0s_vector, EoR_vector)

    # cons
    cons = np.sum(Xfs_vector) - 1

    return cons

constraint6(x, params)

Sets a bound for sum of all mole fraction (x[i] = 1) representing f(x)=1

Parameters

x : list loop values params : list parameters

Returns

cons : float constraint

Notes

1. sum of all mole fraction should be 1

Source code in PyReactLab/docs/optim.py

def constraint6(self, x, params):
    '''
    Sets a bound for sum of all mole fraction (x[i] = 1) representing f(x)=1

    Parameters
    ------------
    x : list
    loop values
    params : list
    parameters

    Returns
    -------
    cons : float
    constraint

    Notes
    -----
    1. sum of all mole fraction should be 1
    '''
    # print(f"5: {x}")
    N0s, comp_list, component_dict = params

    # build EoR
    comp_value_list, comp_value_matrix = self.build_EoR(x)

    # extent sun [mol]
    EoR_vector = np.sum(comp_value_matrix, axis=0)

    # unpack N0s
    N0s_list, N0s_vector, N0f = self.unpack_X(N0s)

    # build final X
    Xfs, Xfs_vector, Nf, Nfs_vector = self.build_final_X(
        N0s_vector, EoR_vector)

    # cons
    cons = np.sum(Xfs_vector) - 1

    return cons

constraints_collection_1(initial_mole)

Generate a collection of constraints for optimization.

Parameters

initial_mole : Dict[str, float | int] Initial mole dictionary {key: value}, such as {CO2: 0, H2: 1, CO: 2, H2O: 3, CH3OH: 4}

Source code in PyReactLab/docs/optim.py

def constraints_collection_1(self, initial_mole: Dict[str, float | int]):
    """
    Generate a collection of constraints for optimization.

    Parameters
    ----------
    initial_mole : Dict[str, float | int]
        Initial mole dictionary {key: value}, such as {CO2: 0, H2: 1, CO: 2, H2O: 3, CH3OH: 4}
    """
    try:
        # NOTE: define constraint
        cons = []

        # inequality
        # looping through components
        for key, value in self.component_dict.items():
            # append constraint
            cons.append({
                'type': 'ineq',
                'fun': self.constraint1,
                'args': (
                    (initial_mole, self.comp_list, self.component_dict, value),)
            })
            # append constraint
            cons.append({
                'type': 'ineq',
                'fun': self.constraint2,
                'args': (
                    (initial_mole, self.comp_list, self.component_dict, value),)
            })
            # append constraint
            cons.append({
                'type': 'ineq',
                'fun': self.constraint3,
                'args': (
                    (initial_mole, self.comp_list, self.component_dict, value),)
            })

            # check consumed
            if key in self.overall_reaction_analysis['consumed']:
                # append constraint
                cons.append({
                    'type': 'ineq',
                    'fun': self.constraint4,
                    'args': (
                        (initial_mole, self.comp_list, self.component_dict, value),)
                })

        # append constraint
        cons.append({
            'type': 'ineq',
            'fun': self.constraint5,
                    'args': ((initial_mole, self.comp_list, self.component_dict),)
        })

        # equality
        # constraint 5
        # cons.append({
        # 'type': 'ineq',
        # 'fun': constraint5,
        # 'args':((N0s, comp_list,component_dict),)
        # })

        # return
        return cons
    except Exception as e:
        raise Exception(
            f"Error in generating constraints collection: {str(e)}") from e

constraints_collection_2(nu, n0, include_mole_fraction_constraint=True)

Generate a collection of constraints for optimization.

Parameters

initial_mole : Dict[str, float | int] Initial mole dictionary {key: value}, such as {CO2: 0, H2: 1, CO: 2, H2O: 3, CH3OH: 4}

Source code in PyReactLab/docs/optim.py

def constraints_collection_2(
    self,
    nu,
    n0,
    include_mole_fraction_constraint=True
):
    """
    Generate a collection of constraints for optimization.

    Parameters
    ----------
    initial_mole : Dict[str, float | int]
        Initial mole dictionary {key: value}, such as {CO2: 0, H2: 1, CO: 2, H2O: 3, CH3OH: 4}
    """
    try:
        n_species, n_reactions = nu.shape
        constraints = []

        # 1. Species non-negativity constraints
        for i in range(n_species):
            def make_fn(i):
                return lambda ksi: n0[i] + np.dot(nu[i, :], ksi)
            constraints.append(NonlinearConstraint(make_fn(i), 0, np.inf))

        # 2. Mole fraction sum constraint (optional)
        if include_mole_fraction_constraint:
            def mole_fraction_sum(ksi):
                n = n0 + nu @ ksi
                total = np.sum(n)
                # avoid division by zero if total is very small
                if total <= 1e-12:
                    return 0.0
                y = n / total
                return np.sum(y)

            constraints.append(NonlinearConstraint(
                mole_fraction_sum, 1.0, 1.0))

        return constraints
    except Exception as e:
        raise Exception(
            f"Error in generating constraints collection: {str(e)}") from e

equilibrium_results(initial_mole, EoR)

Calculate the equilibrium results based on the initial mole and extent of reaction.

Parameters

initial_mole : dict Initial mole dictionary {key: value}, such as {CO2: 0, H2: 1, CO: 2, H2O: 3, CH3OH: 4} EoR : list Extent of reaction list [EoR1, EoR2, EoR3, ...]

Returns

Xfs : dict Final mole fraction dictionary {key: value}, such as {CO2: 0.1, H2: 0.1, CO: 0.2, H2O: 0.3, CH3OH: 0.4} Xfs_vector : np.ndarray Final mole fraction vector [Xfs1, Xfs2, Xfs3, ...] Nf : float Final moles of the system Nfs_vector : np.ndarray Final moles vector [Nf1, Nf2, Nf3, ...] Nfs : dict Final moles dictionary {key: value}, such as {CO2: 0, H2: 1, CO: 2, H2O: 3, CH3OH: 4}

Source code in PyReactLab/docs/optim.py

def equilibrium_results(
    self,
    initial_mole: Dict[str, float | int],
    EoR: list
):
    '''
    Calculate the equilibrium results based on the initial mole and extent of reaction.

    Parameters
    ----------
    initial_mole : dict
        Initial mole dictionary {key: value}, such as {CO2: 0, H2: 1, CO: 2, H2O: 3, CH3OH: 4}
    EoR : list
        Extent of reaction list [EoR1, EoR2, EoR3, ...]

    Returns
    -------
    Xfs : dict
        Final mole fraction dictionary {key: value}, such as {CO2: 0.1, H2: 0.1, CO: 0.2, H2O: 0.3, CH3OH: 0.4}
    Xfs_vector : np.ndarray
        Final mole fraction vector [Xfs1, Xfs2, Xfs3, ...]
    Nf : float
        Final moles of the system
    Nfs_vector : np.ndarray
        Final moles vector [Nf1, Nf2, Nf3, ...]
    Nfs : dict
        Final moles dictionary {key: value}, such as {CO2: 0, H2: 1, CO: 2, H2O: 3, CH3OH: 4}
    '''
    try:
        # unpack
        N0s_list, N0s_vector, N0f = self.unpack_X(initial_mole)

        # update comp_list with EoR
        # build EoR => item[key] = value * EoR[i]
        _, comp_value_matrix = self.build_EoR(EoR)

        # extent sun [mol]
        EoR_vector = np.sum(comp_value_matrix, axis=0)

        # build final X
        Xfs, Xfs_vector, Nf, Nfs_vector = self.build_final_X(
            N0s_vector, EoR_vector
        )

        # Nfs
        Nfs = {}
        for key, value in self.component_dict.items():
            Nfs[key] = Nfs_vector[value]

        # set float
        Nfs = {k: float(v) for k, v in Nfs.items()}
        Xfs = {k: float(v) for k, v in Xfs.items()}

        return Xfs, Xfs_vector, Nf, Nfs_vector, Nfs
    except Exception as e:
        raise Exception(
            f"Error in processing optimization results: {str(e)}") from e

obj_fn(x, N0s, P, T, equilibrium_constants, method)

Objective function for optimization

Parameters

x : list loop values N0s : dict initial mole dictionary {key: value}, such as {CO2: 0, H2: 1, CO: 2, H2O: 3, CH3OH: 4} parameters P : float pressure [bar] T : float temperature [K] equilibrium_constants : dict reaction equilibrium constant dictionary calculated at T as: {key: value}, such as {R1: 0.1, R2: 0.2, ...} method : str optimization method.

Returns

obj : float | list objective function, if method is 'minimize' return float, if method is 'least_squares' return list

Source code in PyReactLab/docs/optim.py

def obj_fn(
    self,
    x,
    N0s: Dict[str, float],
    P: float,
    T: float,
    equilibrium_constants: Dict[str, float],
    method: Literal['minimize', 'least_squares'],
):
    '''
    Objective function for optimization

    Parameters
    ----------
    x : list
        loop values
    N0s : dict
        initial mole dictionary {key: value}, such as {CO2: 0, H2: 1, CO: 2, H2O: 3, CH3OH: 4}
        parameters
    P : float
        pressure [bar]
    T : float
        temperature [K]
    equilibrium_constants : dict
        reaction equilibrium constant dictionary calculated at T as: {key: value}, such as {R1: 0.1, R2: 0.2, ...}
    method : str
        optimization method.

    Returns
    -------
    obj : float | list
        objective function, if method is 'minimize' return float, if method is 'least_squares' return list
    '''
    try:
        # NOTE: default values
        # pressure [bar]
        # self.P_Ref_bar

        # NOTE: extent of reaction
        EoR = x

        # extent of reaction
        EoR = x
        # define threshold
        for i, item in enumerate(EoR):
            EoR[i] = max(self.threshold, item)

        # SECTION: mole balance
        # NOTE: initial mole
        # unpack N0s
        N0s_list, N0s_vector, N0f = self.unpack_X(N0s)

        # NOTE: update comp_list with EoR
        # build EoR => item[key] = value * EoR[i]
        _, comp_value_matrix = self.build_EoR(EoR)

        # NOTE: extent sun [mol]
        EoR_vector = np.sum(comp_value_matrix, axis=0)

        # NOTE: build final X
        Xfs, Xfs_vector, Nf, Nfs_vector = self.build_final_X(
            N0s_vector, EoR_vector)

        # SECTION: component state analysis

        # SECTION: equilibrium equation
        # reaction term
        reaction_terms = self.reaction_equilibrium_equation(
            Xfs=Xfs,
            P=P,
            T=T,
            equilibrium_constants=equilibrium_constants,
        )

        # SECTION: build objective functions
        if method == 'minimize':
            obj = 0
            for i, reaction in enumerate(self.reaction_analysis):
                # NOTE: method 1
                # obj = abs(reaction_terms[reaction])
                # NOTE: method 2
                # obj += reaction_terms[reaction]**2
                # NOTE: set
                obj += reaction_terms[reaction]**2

            # NOTE: define objective function
            # penalty obj
            # obj += 1e-3

            # sqrt
            # obj = sqrt(obj)
        elif method == 'least_squares':
            # NOTE: set
            obj = []
            for i, reaction in enumerate(self.reaction_analysis):
                # NOTE: method 1
                # obj.append(reaction_terms[reaction])
                # NOTE: method 2
                # obj.append(reaction_terms[reaction]**2)
                # NOTE: set
                obj.append(reaction_terms[reaction])

            # NOTE: mole fraction constraints
            Xfs_sum = np.sum(Xfs_vector)
            obj.append(Xfs_sum - 1)
        else:
            raise ValueError(
                f"Invalid optimization method: {method}. Must be 'minimize' or 'least_squares'.")

        return obj
    except Exception as e:
        raise Exception(
            f"Error in objective function: {str(e)}") from e

opt_run(initial_mole, initial_mole_fraction, pressure, temperature, equilibrium_constants, reaction_numbers, method='minimize', **kwargs)

Start optimization process

Parameters

initial_mole : dict Initial mole dictionary {key: value}, such as {CO2: 0, H2: 1, CO: 2, H2O: 3, CH3OH: 4} initial_mole_fraction : dict Initial mole fraction dictionary {key: value}, such as {CO2: 0, H2: 1, CO: 2, H2O: 3, CH3OH: 4} pressure : float Pressure [bar] temperature : float Temperature [K] equilibrium_constants : dict Reaction equilibrium constant dictionary calculated at T as: {key: value}, such as {R1: 0.1, R2: 0.2, ...} reaction_numbers : int Number of reactions method : str, optional Optimization method, by default 'minimize'. kwargs : dict Additional parameters for optimization. - least_square_algorithm: 'trf', 'dogbox', 'lm' - minimize_algorithm: 'SLSQP', 'trust-constr', 'COBYLA' - scale: Scale factor for the upper bound. - bound_scale: Scale factor for the upper bound. - fallback: Fallback value for upper bound if no reactants are present.

Returns

opt_res : OptimizeResult Optimization result object containing the optimized values and other information.

Source code in PyReactLab/docs/optim.py

def opt_run(
    self,
    initial_mole: Dict[str, float | int],
    initial_mole_fraction: Dict[str, float | int],
    pressure: float,
    temperature: float,
    equilibrium_constants: Dict[str, float],
    reaction_numbers: int,
    method: Literal['minimize', 'least_squares'] = 'minimize',
    **kwargs
):
    """
    Start optimization process

    Parameters
    ----------
    initial_mole : dict
        Initial mole dictionary {key: value}, such as {CO2: 0, H2: 1, CO: 2, H2O: 3, CH3OH: 4}
    initial_mole_fraction : dict
        Initial mole fraction dictionary {key: value}, such as {CO2: 0, H2: 1, CO: 2, H2O: 3, CH3OH: 4}
    pressure : float
        Pressure [bar]
    temperature : float
        Temperature [K]
    equilibrium_constants : dict
        Reaction equilibrium constant dictionary calculated at T as: {key: value}, such as {R1: 0.1, R2: 0.2, ...}
    reaction_numbers : int
        Number of reactions
    method : str, optional
        Optimization method, by default 'minimize'.
    kwargs : dict
        Additional parameters for optimization.
        - least_square_algorithm: 'trf', 'dogbox', 'lm'
        - minimize_algorithm: 'SLSQP', 'trust-constr', 'COBYLA'
        - scale: Scale factor for the upper bound.
        - bound_scale: Scale factor for the upper bound.
        - fallback: Fallback value for upper bound if no reactants are present.

    Returns
    -------
    opt_res : OptimizeResult
        Optimization result object containing the optimized values and other information.
    """
    try:
        # NOTE: default values
        # least_square_algorithm
        least_square_algorithm = kwargs.get(
            'least_square_algorithm', 'trf')
        # minimize_algorithm
        minimize_algorithm = kwargs.get('minimize_algorithm', 'SLSQP')

        # ? check if the algorithm is valid
        if minimize_algorithm not in ['SLSQP', 'trust-constr', 'COBYLA']:
            raise ValueError(
                f"Invalid minimize algorithm: {minimize_algorithm}. Must be 'SLSQP', 'trust-constr', or 'COBYLA'.")

        # ? check if the algorithm is valid
        if least_square_algorithm not in ['trf', 'dogbox', 'lm']:
            raise ValueError(
                f"Invalid least square algorithm: {least_square_algorithm}. Must be 'trf', 'dogbox', or 'lm'.")

        # NOTE: set
        # initial guess for extent of reaction
        # default value
        EOR0_val = kwargs.get('EOR0_val', 0.5)

        # set
        EOR0 = np.random.uniform(0, EOR0_val, reaction_numbers)

        # NOTE: bounds
        bound0 = (0, 20)
        bounds = []
        for i in range(len(EOR0)):
            bounds.append(bound0)

        # EOR02
        # initial mole
        initial_mole_ = np.array(list(initial_mole.values()))

        # extent of reaction (EoR)
        EOR0_Bounds, EOR0_bounds, EoR_initial = self.compute_bounds(
            nu=self.stoichiometric_coeff,
            n0=initial_mole_,
            **kwargs,
        )

        # NOTE: constraints
        # constraints 1
        cons1 = self.constraints_collection_1(initial_mole)
        # constraints 2
        cons2 = self.constraints_collection_2(
            nu=self.stoichiometric_coeff,
            n0=initial_mole_,
            include_mole_fraction_constraint=True
        )

        # set constraints
        cons = []

        # SECTION: optimize
        if method == 'minimize':
            opt_res = optimize.minimize(
                fun=self.obj_fn,
                x0=EOR0,
                args=(
                    initial_mole,
                    pressure,
                    temperature,
                    equilibrium_constants,
                    method
                ),
                method=minimize_algorithm,
                bounds=EOR0_bounds,
                constraints=cons1,
                options={
                    'disp': False,
                    'ftol': 1e-12,
                    'maxiter': 1000
                }
            )
        elif method == 'least_squares':
            opt_res = optimize.least_squares(
                fun=self.obj_fn,
                x0=EOR0,
                args=(
                    initial_mole,
                    pressure,
                    temperature,
                    equilibrium_constants,
                    method
                ),
                bounds=EOR0_Bounds,
                method=least_square_algorithm,
            )
        else:
            raise ValueError(
                f"Invalid optimization method: {method}. Must be 'minimize' or 'least_squares'.")

        # save
        return opt_res
    except Exception as e:
        raise Exception(
            f"Error in the optimization process: {str(e)}") from e

reaction_equilibrium_equation(Xfs, P, T, equilibrium_constants, phase='gas', **kwargs)

Generate reaction equilibrium equations for gas and liquid phases.

Parameters

Xfs : dict Mole fraction dictionary {key: value}, such as {CO2: 0, H2: 1, CO: 2, H2O: 3, CH3OH: 4} P : float Pressure [Pa] T : float Temperature [K] equilibrium_constants : dict Reaction equilibrium constant dictionary calculated at T as: such as {R1: dict, R2: dict, ...}, the dict value consists of value, symbol, unit, temperature, reaction, method phase : str, optional Phase of calculation, by default "gas". kwargs : dict Additional parameters for non-ideal calculations.

Notes

The equilibrium constant at temperature T is calculated from the standard Gibbs free energy change:

K_T = exp(-ΔG_T⁰ / RT) = ∏ [ (f̂_i / f_i⁰) ^ ν_i ]

Where: - K_T : Equilibrium constant at temperature T - ΔG_T⁰ : Standard Gibbs free energy change at temperature T (J/mol) - R : Universal gas constant (8.314 J/mol·K) - T : Temperature in Kelvin - f̂_i : The mixture fugacity of component i. - f_i⁰ : The fugacity of pure component i in its reference state at temperature T. - ν_i : Stoichiometric coefficient of component i (positive for products, negative for reactants)

This expression connects thermodynamic data to fugacity-based equilibrium conditions.

Gas Phase Calculation: - For gaseous compounds, the reference is the ideal gas at 1 bar, then f_i⁰ = P_ref = 1 bar. - f̂_i for non-ideal gases is calculated as f̂_i = φ_i * P * Y_i. - φ_i is the fugacity coefficient in the mixture.

Liquid Phase Calculation: - f_i⁰ is standard-state fugacity (e.g., pure liquid at 1 bar) - f̂_i for liquid mixtures is defined as f̂_i = γ_i * X_i * f_i(pure,T). - f_i(pure,T) is the fugacity of pure component i at temperature T. - γ_i is the activity coefficient of component i in the liquid phase. - X_i is the mole fraction of component i in the liquid phase.

Assumptions: - Regular liquid mixtures: Lewis-Randall - Ideal liquid solution: Raoult - Dilute or electrolyte systems: Henry - Poynting correction for liquid mixtures at high pressures. - General theoretical frameworks: Fugacity-based (chemical potential)

Source code in PyReactLab/docs/optim.py

def reaction_equilibrium_equation(
    self,
    Xfs: Dict[str, float],
    P: float,
    T: float,
    equilibrium_constants: Dict[str, Dict[str, Any]],
    phase: Literal[
        "liquid", "gas"
    ] = "gas",
    **kwargs
):
    """
    Generate reaction equilibrium equations for gas and liquid phases.

    Parameters
    ----------
    Xfs : dict
        Mole fraction dictionary {key: value}, such as {CO2: 0, H2: 1, CO: 2, H2O: 3, CH3OH: 4}
    P : float
        Pressure [Pa]
    T : float
        Temperature [K]
    equilibrium_constants : dict
        Reaction equilibrium constant dictionary calculated at T as: such as {R1: dict, R2: dict, ...},
        the dict value consists of `value`, `symbol`, `unit`, `temperature`, `reaction`, `method`
    phase : str, optional
        Phase of calculation, by default "gas".
    kwargs : dict
        Additional parameters for non-ideal calculations.

    Notes
    -----
    The equilibrium constant at temperature T is calculated from the standard Gibbs free energy change:

    K_T = exp(-ΔG_T⁰ / RT) = ∏ [ (f̂_i / f_i⁰) ^ ν_i ]

    Where:
    - K_T    : Equilibrium constant at temperature T
    - ΔG_T⁰  : Standard Gibbs free energy change at temperature T (J/mol)
    - R      : Universal gas constant (8.314 J/mol·K)
    - T      : Temperature in Kelvin
    - f̂_i    : The mixture fugacity of component i.
    - f_i⁰   : The fugacity of pure component i in its reference state at temperature T.
    - ν_i    : Stoichiometric coefficient of component i (positive for products, negative for reactants)

    This expression connects thermodynamic data to fugacity-based equilibrium conditions.

    Gas Phase Calculation:
    - For gaseous compounds, the reference is the ideal gas at 1 bar, then f_i⁰ = P_ref = 1 bar.
    - f̂_i for non-ideal gases is calculated as f̂_i = φ_i * P * Y_i.
    - φ_i is the fugacity coefficient in the mixture.

    Liquid Phase Calculation:
    - f_i⁰ is standard-state fugacity (e.g., pure liquid at 1 bar)
    - f̂_i for liquid mixtures is defined as f̂_i = γ_i * X_i * f_i(pure,T).
    - f_i(pure,T) is the fugacity of pure component i at temperature T.
    - γ_i is the activity coefficient of component i in the liquid phase.
    - X_i is the mole fraction of component i in the liquid phase.

    Assumptions:
    - Regular liquid mixtures: Lewis-Randall
    - Ideal liquid solution: Raoult
    - Dilute or electrolyte systems: Henry
    - Poynting correction for liquid mixtures at high pressures.
    - General theoretical frameworks: Fugacity-based (chemical potential)
    """
    try:
        # NOTE: default values
        # pressure [bar]
        # self.P_Ref_bar
        # temperature [K]
        # self.T_Ref_K

        # SECTION: fugacity coefficient
        # fugacity coefficient
        fugacity_coeff = {}

        # check gas mixture
        if self.gas_mixture.lower() == "non-ideal":
            # NOTE: model input
            model_input = {
                "feed-specification": Xfs,
                "pressure": [P, 'bar'],
                "temperature": [T, 'K'],
            }

            # NOTE: eos model
            eos_model = self.eos_model

            # Validate the model name
            if eos_model not in EOS_MODELS:
                raise ValueError(
                    f"Invalid EOS model: {eos_model}. Must be {EOS_MODELS}.")

            # NOTE: calculate fugacity
            res_ = self._cal_fugacity_coefficient_gaseous_mixture(
                model_name=eos_model,  # type: ignore
                model_input=model_input)
            # update
            fugacity_coeff = {**res_}

        elif self.gas_mixture.lower() == "ideal":
            # set
            for i, key in enumerate(self.component_dict.keys()):
                fugacity_coeff[key] = 1
        else:
            raise ValueError(
                "Invalid gas mixture mode. Must be 'ideal' or 'non-ideal'.")

        # SECTION: activity coefficient
        # activity coefficient
        activity_coeff = {}

        # check liquid mixture
        if self.solution.lower() == "non-ideal":

            # prepare model input
            # Ensure activity_inputs is a dictionary before unpacking
            activity_inputs = self.activity_inputs or {}
            # set
            model_inputs = {
                "mole_fraction": Xfs, **activity_inputs
            }

            # NOTE: check model name
            if self.activity_model == 'NRTL':
                # exec
                activity_coeff = self._cal_activity_coefficient_solution(
                    model_name='NRTL', model_input=model_inputs)
            elif self.activity_model == 'UNIQUAC':
                # exec
                activity_coeff = self._cal_activity_coefficient_solution(
                    model_name='UNIQUAC', model_input=model_inputs)
            else:
                raise ValueError(
                    f"Invalid activity model: {self.activity_model}. Must be {ACTIVITY_MODELS}.")
        elif self.solution.lower() == "ideal":
            # set
            for i, key in enumerate(self.component_dict.keys()):
                activity_coeff[key] = 1
        else:
            raise ValueError(
                "Invalid liquid mixture mode. Must be 'ideal' or 'non-ideal'.")

        # SECTION: equilibrium equation
        # reaction term
        reaction_terms = {}

        # NOTE: loop over reactions
        for i, reaction in enumerate(self.reaction_analysis):
            # denominator
            denominator = 1
            # numerator
            numerator = 1

            # state count
            state_count_ = self.reaction_analysis[reaction]['state_count']

            # SECTION: loop over reactants
            for item in self.reaction_analysis[reaction]['reactants']:
                # NOTE: item info
                molecule_ = item['molecule']
                molecule_state_ = item['molecule_state']
                coefficient_ = item['coefficient']
                state_ = item['state']

                # final mole fraction
                Xfs_ = Xfs[molecule_state_]

                # NOTE: cal
                # check phase
                if state_ == "g":
                    # set
                    term_ = Xfs_ * (P / self.P_Ref_bar) * \
                        fugacity_coeff[molecule_state_]
                elif state_ == "l":
                    # check solution
                    if self.solution.lower() == "non-ideal":
                        # Lewis-Randall/Raoult
                        term_ = Xfs_ * activity_coeff[molecule_state_]
                    elif self.solution.lower() == "ideal":
                        # check
                        if state_count_['l'] == 1:
                            # pure liquid
                            term_ = 1
                    else:
                        raise ValueError(
                            "Invalid liquid mixture mode. Must be 'ideal' or 'non-ideal'.")

                elif state_ == "s":
                    # set
                    # solid always has the activity of 1
                    term_ = 1
                else:
                    raise ValueError(
                        "Invalid phase. Must be 'gas' or 'liquid'.")

                # update denominator
                denominator *= (term_)**coefficient_

            # SECTION: loop over products
            for item in self.reaction_analysis[reaction]['products']:
                # NOTE: item info
                molecule_ = item['molecule']
                molecule_state_ = item['molecule_state']
                coefficient_ = item['coefficient']
                state_ = item['state']

                # final mole fraction
                Xfs_ = Xfs[molecule_state_]

                # NOTE: cal
                # check phase
                if state_ == "g":
                    # set
                    term_ = Xfs_ * (P / self.P_Ref_bar) * \
                        fugacity_coeff[molecule_state_]
                elif state_ == "l":
                    # check solution
                    if self.solution.lower() == "non-ideal":
                        # Lewis-Randall/Raoult
                        term_ = Xfs_ * activity_coeff[molecule_state_]
                    elif self.solution.lower() == "ideal":
                        # check
                        if state_count_['l'] == 1:
                            # pure liquid
                            term_ = 1
                    else:
                        raise ValueError(
                            "Invalid liquid mixture mode. Must be 'ideal' or 'non-ideal'.")
                elif state_ == "s":
                    # set
                    # solid always has the activity of 1
                    term_ = 1
                else:
                    raise ValueError(
                        "Invalid phase. Must be 'gas' or 'liquid'.")

                # update numerator
                numerator *= (term_)**coefficient_

            # NOTE: update reaction term
            # equilibrium data
            Keq_ = equilibrium_constants[reaction]['value']

            # ? method 1
            # reaction_terms[reaction] = (numerator/denominator) - Keq_
            # ? method 2
            # reaction_terms[reaction] = numerator - denominator*Keq_
            # ? method 3
            Q_i = numerator / denominator
            Q_i_safe = np.clip(Q_i, 1e-12, None)
            reaction_terms[reaction] = log(Q_i_safe) - log(Keq_)

        # return
        return reaction_terms
    except Exception as e:
        raise Exception(
            f"Error in building the reaction equilibrium equation for {phase}: {str(e)}") from e

unpack_X(mol_data_pack)

convert X to dict {key: value} and list [value]

Parameters

mol_data_pack : dict component mole dictionary {key: value}, such as {CO2: 0, H2: 1, CO: 2, H2O: 3, CH3OH: 4}

Returns

N0s_list : list The list of N0s N0s_vector : numpy.ndarray The vector of N0s Nf : float The total number of moles

Source code in PyReactLab/docs/optim.py

def unpack_X(self, mol_data_pack: Dict[str, float | int]):
    '''
    convert X to dict {key: value} and list [value]

    Parameters
    ----------
    mol_data_pack : dict
        component mole dictionary {key: value}, such as {CO2: 0, H2: 1, CO: 2, H2O: 3, CH3OH: 4}

    Returns
    -------
    N0s_list : list
        The list of N0s
    N0s_vector : numpy.ndarray
        The vector of N0s
    Nf : float
        The total number of moles
    '''
    try:
        # size
        size = len(mol_data_pack)

        # lists
        N0s_list = []
        N0s_vector = np.zeros(size)

        # looping through
        for i, (com_key, com_value) in enumerate(self.component_dict.items()):
            N0s_list.append(mol_data_pack[com_key])
            N0s_vector[i] = mol_data_pack[com_key]

        # final
        Nf = np.sum(N0s_vector)

        return N0s_list, N0s_vector, Nf
    except Exception as e:
        raise Exception(
            f"Error in ReactionOptimizer.unpack_X(): {str(e)}") from e