"""This module for defining chemical reaction objects was originally sourced from
pymatgen and streamlined for the reaction-network code.
"""
from __future__ import annotations
import re
from copy import deepcopy
from functools import cached_property
from itertools import chain, combinations
from typing import TYPE_CHECKING
import numpy as np
from monty.fractions import gcd_float
from rxn_network.core import Composition
from rxn_network.data import COMMON_GASES
from rxn_network.reactions.base import Reaction
if TYPE_CHECKING:
from collections.abc import Iterable
from pymatgen.core.periodic_table import Element
TOLERANCE = 1e-6 # Tolerance for determining if a particular component fraction is > 0.
[docs]
class BasicReaction(Reaction):
"""An object representing a basic chemical reaction: compositions and their
coefficients.
"""
def __init__(
self,
compositions: Iterable[Composition],
coefficients: Iterable[float] | np.ndarray,
balanced: bool | None = None,
data: dict | None = None,
lowest_num_errors: float = 0,
):
"""A BasicReaction object is defined by a list of compositions and their
corresponding coefficients, where a negative coefficient refers to a
reactant, and a positive coefficient refers to a product.
Args:
compositions: List of composition objects (pymatgen).
coefficients: List of coefficients, where negative coeff distinguishes a
reactant.
balanced: Whether the reaction is stoichiometricaly balanced or not
(see construction via balance() method).
data: Optional corresponding data in dictionary format; often used to store
various calculated parameters.
lowest_num_errors: the minimum number of errors reported by the reaction
balancing algorithm (see the balance() method). A number of errors
>= 1 means that the reaction may be different than intended (some
phases may be shuffled or removed entirely).
"""
self._compositions = [Composition(c) for c in compositions]
self._coefficients = np.array(coefficients)
self.reactant_coeffs = {comp: coeff for comp, coeff in zip(self._compositions, self._coefficients) if coeff < 0}
self.product_coeffs = {comp: coeff for comp, coeff in zip(self._compositions, self._coefficients) if coeff > 0}
if balanced is not None:
self.balanced = balanced
else:
sum_reactants = sum([k * abs(v) for k, v in self.reactant_coeffs.items()], Composition({}))
sum_products = sum([k * abs(v) for k, v in self.product_coeffs.items()], Composition({}))
if not sum_reactants.almost_equals(sum_products, rtol=0, atol=TOLERANCE):
self.balanced = False
else:
self.balanced = True
self.data = data if data else {}
self.lowest_num_errors = lowest_num_errors
[docs]
@classmethod
def balance(
cls,
reactants: list[Composition],
products: list[Composition],
data: dict | None = None,
) -> BasicReaction:
"""Reactants and products to be specified as list of
pymatgen.core.Composition. e.g., [comp1, comp2].
Args:
reactants: List of reactants.
products: List of products.
data: Optional dictionary containing extra data about the reaction.
"""
compositions = list(reactants + products)
coeffs, lowest_num_errors, num_constraints = cls._balance_coeffs(reactants, products)
if not data:
data = {}
data["num_constraints"] = num_constraints
balanced = True
if coeffs is None or lowest_num_errors == np.inf:
balanced = False
coeffs = np.zeros(len(compositions))
return cls(
compositions=compositions,
coefficients=coeffs,
balanced=balanced,
data=data,
lowest_num_errors=lowest_num_errors,
)
[docs]
def normalize_to(self, comp: Composition, factor: float = 1) -> BasicReaction:
"""Normalizes the reaction to one of the compositions via the provided factor.
By default, normalizes such that the composition given has a coefficient of
1.
Args:
comp: Composition object to normalize to
factor: factor to normalize to. Defaults to 1.
"""
all_comp = self.compositions
coeffs = self.coefficients.copy()
scale_factor = abs(1 / coeffs[self.compositions.index(comp)] * factor)
coeffs *= scale_factor
return BasicReaction(all_comp, coeffs)
[docs]
def normalize_to_element(self, element: Element, factor: float = 1) -> BasicReaction:
"""Normalizes the reaction to one of the elements.
By default, normalizes such that the amount of the element is 1.
Another factor can be specified.
Args:
element (Element/Species): Element to normalize to.
factor (float): Factor to normalize to. Defaults to 1.
"""
all_comp = self.compositions
coeffs = self.coefficients.copy()
current_el_amount = sum(all_comp[i][element] * abs(coeffs[i]) for i in range(len(all_comp))) / 2
scale_factor = factor / current_el_amount
coeffs *= scale_factor
return BasicReaction(all_comp, coeffs)
[docs]
def get_el_amount(self, element: Element) -> float:
"""Returns the amount of the element in the reaction.
Args:
element (Element/Species): Element in the reaction
Returns:
Amount of that element in the reaction.
"""
return sum(self.compositions[i][element] * abs(self.coefficients[i]) for i in range(len(self.compositions))) / 2
[docs]
def get_coeff(self, comp: Composition):
"""Returns coefficient for a particular composition."""
return self.coefficients[self.compositions.index(comp)]
[docs]
def normalized_repr_and_factor(self):
"""Normalized representation for a reaction
For example, ``4 Li + 2 O -> 2Li2O`` becomes ``2 Li + O -> Li2O``.
"""
return self._str_from_comp(self.coefficients, self.compositions, reduce=True)
[docs]
def copy(self) -> BasicReaction:
"""Returns a copy of the BasicReaction object."""
return BasicReaction(
compositions=self.compositions,
coefficients=self.coefficients.copy(),
balanced=self.balanced,
data=self.data,
lowest_num_errors=self.lowest_num_errors,
)
[docs]
def reverse(self) -> BasicReaction:
"""Returns a copy of the original BasicReaction object where original reactants are
new products, and vice versa.
"""
return BasicReaction(
compositions=self.compositions,
coefficients=-1 * self.coefficients.copy(),
balanced=self.balanced,
data=self.data,
lowest_num_errors=self.lowest_num_errors,
)
[docs]
def is_separable(self, target: str | Composition) -> bool:
"""Checks if the reaction forms byproducts which are separable from the target
composition.
Separable byproducts are those that are common gases (e.g., CO2), or other phases that do not contain any of the
elements in the target phase.
Args:
target: Composition of target; elements in this phase will be used to
determine whether byproducts only contain added elements.
Returns:
True if reaction is separable from target, False otherwise.
"""
target = Composition(target)
identified_targets = [c for c in self.compositions if c.reduced_composition == target.reduced_composition]
if len(identified_targets) == 0:
raise ValueError(f"Target composition {target} not in reaction {self}")
added_elems = set(self.elements) - set(target.elements)
products = set(deepcopy(self.products))
for t in identified_targets:
products.remove(t)
are_separable = []
for comp in products:
if comp in COMMON_GASES or added_elems.issuperset(comp.elements):
are_separable.append(True)
else:
are_separable.append(False)
return all(are_separable)
@cached_property
def reactant_atomic_fractions(self) -> dict:
"""Returns the atomic mixing ratio of reactants in the reaction."""
if not self.balanced:
raise ValueError("Reaction is not balanced")
return {
c.reduced_composition: -coeff * c.num_atoms / self.num_atoms for c, coeff in self.reactant_coeffs.items()
}
@cached_property
def product_atomic_fractions(self) -> dict:
"""Returns the atomic mixing ratio of reactants in the reaction."""
if not self.balanced:
raise ValueError("Reaction is not balanced")
return {c.reduced_composition: coeff * c.num_atoms / self.num_atoms for c, coeff in self.product_coeffs.items()}
@cached_property
def reactant_molar_fractions(self) -> dict:
"""Returns the molar mixing ratio of reactants in the reaction."""
if not self.balanced:
raise ValueError("Reaction is not balanced")
total = sum(self.reactant_coeffs.values())
return {c: coeff / total for c, coeff in self.reactant_coeffs.items()}
@cached_property
def product_molar_fractions(self) -> dict:
"""Returns the molar mixing ratio of products in the reaction."""
if not self.balanced:
raise ValueError("Reaction is not balanced")
total = sum(self.product_coeffs.values())
return {c: coeff / total for c, coeff in self.product_coeffs.items()}
[docs]
@classmethod
def from_string(cls, rxn_string: str) -> BasicReaction:
"""Generates a balanced reaction from a string. The reaction must
already be balanced.
Args:
rxn_string:
The reaction string. For example, "4 Li + O2-> 2 Li2O"
Returns:
BalancedReaction
"""
rct_str, prod_str = rxn_string.split("->")
def get_comp_amt(comp_str):
return {
Composition(m.group(2)): float(m.group(1) or 1)
for m in re.finditer(r"([\d\.]*(?:[eE]-?[\d\.]+)?)\s*([A-Z][\w\.\(\)]*)", comp_str)
}
reactant_coeffs = get_comp_amt(rct_str)
product_coeffs = get_comp_amt(prod_str)
return cls._from_coeff_dicts(reactant_coeffs, product_coeffs)
@property
def reactants(self) -> list[Composition]:
"""List of reactants for this reaction."""
return list(self.reactant_coeffs.keys())
@property
def products(self) -> list[Composition]:
"""List of products for this reaction."""
return list(self.product_coeffs.keys())
@property
def compositions(self) -> list[Composition]:
"""List of composition objects for this reaction."""
return self._compositions
@property
def coefficients(self) -> np.ndarray:
"""Array of reaction coefficients."""
return self._coefficients
@cached_property
def num_atoms(self) -> float:
"""Total number of atoms in this reaction."""
return sum(coeff * sum(c[el] for el in self.elements) for c, coeff in self.product_coeffs.items())
@cached_property
def energy(self) -> float:
"""The energy of this reaction."""
raise ValueError("No energy for a basic reaction!")
@cached_property
def energy_per_atom(self) -> float:
"""The energy per atom of this reaction."""
raise ValueError("No energy per atom for a basic reaction!")
@cached_property
def is_identity(self) -> bool:
"""Returns True if the reaction has identical reactants and products."""
return self._get_is_identity()
def _get_is_identity(self):
"""Returns True if the reaction has identical reactants and products."""
if set(self.reactants) != set(self.products):
return False
if self.balanced is False: # if not balanced, can not check coefficients
return True
return all(np.isclose(self.reactant_coeffs[c] * -1, self.product_coeffs[c]) for c in self.reactant_coeffs)
@cached_property
def chemical_system(self) -> str:
"""Returns the chemical system as string in the form of A-B-C-..."""
return "-".join(sorted([str(el) for el in self.elements]))
@property
def normalized_repr(self) -> str:
"""A normalized representation of the reaction. All factors are converted
to lowest common factors.
"""
return self.normalized_repr_and_factor()[0]
@classmethod
def _balance_coeffs(
cls, reactants: list[Composition], products: list[Composition]
) -> tuple[np.ndarray, int | float, int]:
"""Balances the reaction and returns the new coefficient matrix."""
compositions = reactants + products
num_comp = len(compositions)
all_elems = sorted({elem for c in compositions for elem in c.elements})
num_elems = len(all_elems)
comp_matrix = np.array([[c[el] for el in all_elems] for c in compositions]).T
rank = np.linalg.matrix_rank(comp_matrix)
diff = num_comp - rank
num_constraints = diff if diff >= 2 else 1
# an error = a component changing sides or disappearing
lowest_num_errors = np.inf
first_product_idx = len(reactants)
# start with simplest product constraints, work to more complex constraints
product_constraints = chain.from_iterable(
[combinations(range(first_product_idx, num_comp), n_constr) for n_constr in range(num_constraints, 0, -1)]
)
reactant_constraints = chain.from_iterable(
[combinations(range(first_product_idx), n_constr) for n_constr in range(num_constraints, 0, -1)]
)
best_soln = np.zeros(num_comp)
for constraints in chain(product_constraints, reactant_constraints):
n_constr = len(constraints)
comp_and_constraints = np.append(comp_matrix, np.zeros((n_constr, num_comp)), axis=0)
b = np.zeros((num_elems + n_constr, 1))
b[-n_constr:] = 1 if min(constraints) >= first_product_idx else -1
for num, idx in enumerate(constraints):
comp_and_constraints[num_elems + num, idx] = 1
# arbitrarily fix coeff to 1
coeffs = np.matmul(np.linalg.pinv(comp_and_constraints), b)
num_errors = 0
if np.allclose(np.matmul(comp_matrix, coeffs), np.zeros((num_elems, 1))):
expected_signs = np.array([-1] * len(reactants) + [+1] * len(products))
num_errors = np.sum(np.multiply(expected_signs, coeffs.T) < TOLERANCE)
if num_errors == 0:
lowest_num_errors = 0
best_soln = coeffs
break
if num_errors < lowest_num_errors:
lowest_num_errors = num_errors
best_soln = coeffs
return np.squeeze(best_soln), lowest_num_errors, num_constraints
@staticmethod
def _from_coeff_dicts(reactant_coeffs, product_coeffs) -> BasicReaction:
reactant_comps, r_coefs = zip(*[(comp, -1 * coeff) for comp, coeff in reactant_coeffs.items()])
product_comps, p_coefs = zip(*list(product_coeffs.items()))
return BasicReaction(reactant_comps + product_comps, r_coefs + p_coefs)
@staticmethod
def _str_from_formulas(coeffs, formulas, tol=TOLERANCE) -> str:
reactant_str = []
product_str = []
for amt, formula in zip(coeffs, formulas):
if abs(amt + 1) < tol:
reactant_str.append(formula)
elif abs(amt - 1) < tol:
product_str.append(formula)
elif amt < -tol:
reactant_str.append(f"{-amt:.4g} {formula}")
elif amt > tol:
product_str.append(f"{amt:.4g} {formula}")
return " + ".join(reactant_str) + " -> " + " + ".join(product_str)
@classmethod
def _str_from_comp(cls, coeffs, compositions, reduce=False) -> tuple[str, float]:
r_coeffs = np.zeros(len(coeffs))
r_formulas = []
for i, (amt, comp) in enumerate(zip(coeffs, compositions)):
formula, factor = comp.get_reduced_formula_and_factor()
r_coeffs[i] = amt * factor
r_formulas.append(formula)
if reduce:
factor = 1 / gcd_float(np.abs(r_coeffs))
r_coeffs *= factor
else:
factor = 1
return cls._str_from_formulas(r_coeffs, r_formulas), factor
def __eq__(self, other) -> bool:
if self is other:
return True
if not self.chemical_system == other.chemical_system:
return False
if not len(self.products) == len(other.products):
return False
if not len(self.reactants) == len(other.reactants):
return False
if not np.allclose(sorted(self.coefficients), sorted(other.coefficients)):
return False
if not set(self.reactants) == set(other.reactants):
return False
if not set(self.products) == set(other.products):
return False
return True
def __hash__(self):
return hash(
(self.chemical_system, tuple(sorted(self.coefficients)))
) # not checking here for reactions that are multiples (too expensive)
def __str__(self) -> str:
return self._str_from_comp(self.coefficients, self.compositions)[0]
def __repr__(self) -> str:
return self.__str__()