ChemPy

ChemPy is a Python library designed for researchers and engineers to solve chemistry-related problems through modeling. It provides core abstractions for representing substances, reactions, and reaction systems, and supports both numerical and symbolic computation.

Key Capabilities:

• Chemical Representation: Classes for substances (with composition, weight, charge, and pretty-printing via LaTeX/HTML/Unicode) and reactions (defined by stoichiometry and thermodynamic/kinetic parameters).
• Reaction Systems: Tools to analyze reaction networks (e.g., identifying disjoint sets or connectivity) and express them as:
  • Systems of Ordinary Differential Equations (ODEs) for chemical kinetics.
  • Systems of non-linear equations for chemical equilibria (including multi-phase systems).
• Mathematical Solving:
  • Symbolic Math: Uses SymPy to solve stoichiometry problems analytically, including under-determined systems.
  • Numerical Solving: Integrates with pyodesys for ODEs and pyneqsys for non-linear equations.
  • Optimization: Uses PuLP to express under-determined systems in a canonical form with minimized coefficients.
• Data Support: Supports NumPy arrays, symbolic inputs, and arrays with explicit units for dimensionality checking (e.g., verifying reaction order against rate constants).

ChemPy provides Equilibrium objects to represent chemical equilibria and perform algebraic manipulations like eliminating species or combining reactions.

• Eliminate species: Use Equilibrium.eliminate([eq1, eq2, ...], 'species_name') to find coefficients that cancel out a specific species.
• Combine reactions: Multiply an Equilibrium object by a coefficient and add it to another to create a new combined reaction.
• Cancel species: Use redox.cancel(autoprot) to find a multiplier that allows adding an autoprotolysis reaction (like water) to a redox reaction to simplify the system.

from chempy import Equilibrium
from sympy import symbols

K1, K2 = symbols('K1 K2')
e1 = Equilibrium({'MnO4-': 1, 'H+': 8, 'e-': 5}, {'Mn+2': 1, 'H2O': 4}, K1)
e2 = Equilibrium({'O2': 1, 'H2O': 2, 'e-': 4}, {'OH-': 4}, K2)

# Eliminate 'e-' to combine reactions
coeff = Equilibrium.eliminate([e1, e2], 'e-')
redox = e1*coeff[0] + e2*coeff[1]

To investigate the time dependence of a chemical system (chemical kinetics), use ReactionSystem and get_odesys from chempy.kinetics.ode.

1. Define a ReactionSystem using from_string() with reaction formulas and rate constants.
2. Convert the system to an ODE system using get_odesys(rsys).
3. Integrate the system over a time array using odesys.integrate(tout, c0, ...), where c0 is the initial concentration dictionary.
4. Use the returned result object to plot concentrations over time.

from chempy import ReactionSystem
from chempy.kinetics.ode import get_odesys
import numpy as np

rsys = ReactionSystem.from_string("""2 Fe+2 + H2O2 -> 2 Fe+3 + 2 OH-; 42
    2 Fe+3 + H2O2 -> 2 Fe+2 + O2 + 2 H+; 17
    H+ + OH- -> H2O; 1e10
    H2O -> H+ + OH-; 1e-4""")

odesys, extra = get_odesys(rsys)
tout = np.linspace(0, 23, 100)
c0 = {'Fe+2': 0.05, 'H2O2': 0.1, 'H2O': 1.0, 'H+': 1e-2, 'OH-': 1e-12}

result = odesys.integrate(tout, c0)
result.plot()

The NumSys hierarchy (and its subclasses like NumSysLin, NumSysLog, NumSysSquare, NumSysLinTanh) provides a way to transform the equilibrium equations into a form more suitable for numerical solvers. This is achieved through a transformation pipeline:

1. Pre-processing: Transforms the solver's input variables (e.g., from log-space back to concentration space).
2. Function f(yvec, params): Defines the residual function (the equations that should equal zero) in the transformed space.
3. Post-processing: Transforms the solver's output back into physical concentrations.

Common Transformation Strategies

• NumSysLog: Solves in log-space ($\ln(C)$). This is highly effective for handling concentrations spanning many orders of magnitude and naturally enforces positivity.
• NumSysSquare: Solves for the square of the concentration ($C^2$), which can help with certain algebraic structures.
• NumSysLin: Standard linear transformation approach.
• NumSysLinTanh: Uses a hyperbolic tangent transformation to map the solver's domain to a bounded concentration range $[0, y_{max}]$.

For systems with widely varying concentration scales (stiff problems), you can improve stability using two main approaches:

1. Scaling: Use pyodesys.symbolic.ScaledSys via get_odesys to scale dependencies.
2. Log-space integration: Solve for the logarithm of concentrations to avoid negative values. This is done by creating a custom SymbolicSys using pyodesys.symbolic.symmetricsys with sympy.log and sympy.exp.

Note: Log-space integration can significantly increase the number of function evaluations (computational overhead).

# Log-space integration setup
from pyodesys.symbolic import symmetricsys
import sympy

logexp = (sympy.log, sympy.exp)
LogLogSys = symmetricsys(logexp, logexp, exprs_process_cb=lambda exprs: [sympy.powsimp(expr.expand(), force=True) for expr in exprs])
loglogsys, _ = get_odesys(rsys, SymbolicSys=LogLogSys)

# Integrate using the log-space system
result = loglogsys.integrate(tout, c0, integrator='gsl')