Chempy

Core Principles

Chemical Entities

Represented as strings ('H2O', 'Fe+3') or Substance objects. ChemPy can parse these to determine composition and charge.

Reaction Systems

A collection of Reaction objects. This system can be converted into a mathematical model (ODE) where the rates are defined by mass-action kinetics or custom rate laws.

Units

ChemPy integrates with quantities or pint to ensure dimensional correctness, which is vital for calculating rates (e.g., M⁻¹s⁻¹).

Quick Reference

Installation

pip install chempy

Standard Imports

import numpy as np
import matplotlib.pyplot as plt
from chempy import Reaction, ReactionSystem, Substance
from chempy.util.parsing import formula_to_composition

Basic Pattern - Balancing a Reaction

from chempy import balance_stoichiometry

# 1. Define reactants and products
reac, prod = balance_stoichiometry({'H2', 'O2'}, {'H2O'})

# 2. Output coefficients
print(reac) # {'H2': 2, 'O2': 1}
print(prod) # {'H2O': 2}

Critical Rules

✅ DO

• Use Standard Notation - Follow standard chemical notation for formulas ('H2SO4', 'Fe+3') for reliable parsing.
• Define a ReactionSystem - For any complex network, wrap reactions in a ReactionSystem to manage substances and rates together.
• Check Mass/Charge Balance - Always verify that your reactions are balanced before running kinetic simulations.
• Use Units - Whenever possible, use a units library to avoid errors in time scales or concentration units.
• Specify Rate Laws - Be explicit about whether a reaction is elementary (mass-action) or follows a custom rate law.
• Vectorize Concentrations - Use NumPy arrays when providing initial concentrations to solvers.

❌ DON'T

• Manually Balance Complex Equations - Let balance_stoichiometry handle it, especially for redox reactions.
• Ignore Ionic Strength - In analytical chemistry, remember that activity coefficients change with ionic strength (use chempy.einstein).
• Assume Fast Equilibria - In kinetic models, ensure your rate constants for "instantaneous" steps are high enough but don't cause numerical stiffness.
• Hardcode Molar Masses - Use Substance.from_formula('H2O').mass to ensure precision.

Anti-Patterns (NEVER)

from chempy import balance_stoichiometry

# ❌ BAD: Manual string parsing to find mass
# mass = 1.008 * 2 + 16.00 # Fragile and tedious

# ✅ GOOD: Use Substance properties
from chempy import Substance
water = Substance.from_formula('H2O')
print(water.mass)

# ❌ BAD: Hardcoding ODEs for kinetics
# def dc_dt(c, t): return [-k*c[0]*c[1], ...] # Error-prone

# ✅ GOOD: Generate ODE from ReactionSystem
from chempy.kinetics.ode import get_odesys
rsys = ReactionSystem.from_string("A + B -> C; k")
odesys, extra = get_odesys(rsys)
# This handles all derivatives automatically

Stoichiometry and Formulas

Advanced Balancing (Redox)

from chempy import balance_stoichiometry

# Balancing KMn04 + HCl reaction
reactants = {'KMnO4', 'HCl'}
products = {'KCl', 'MnCl2', 'H2O', 'Cl2'}
reac, prod = balance_stoichiometry(reactants, products)

print(f"Balanced: 2 KMnO4 + 16 HCl -> 2 KCl + 2 MnCl2 + 8 H2O + 5 Cl2")

Substance Properties

from chempy import Substance

# Create substance with metadata
ferric = Substance('Fe+3', name='Iron(III) ion', latex='Fe^{3+}')
print(ferric.composition) # {26: 1} - Atomic number 26
print(ferric.charge)      # 3

Chemical Kinetics (chempy.kinetics)

Simulating Concentration over Time

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

# 1. Define the system: 2A -> B (rate constant k=0.5)
rsys = ReactionSystem.from_string("2 A -> B; 0.5")

# 2. Get the ODE system
odesys, extra = get_odesys(rsys)

# 3. Integrate
tout = np.linspace(0, 10, 50)
c0 = {'A': 1.0, 'B': 0.0}
result = odesys.integrate(tout, c0)

# 4. Plot
plt.plot(result.tout, result.cout)
plt.legend(rsys.substances.keys())

Chemical Equilibria

Calculating Speciation

from chempy.equilibria import EqSystem

# Define equilibria: H2O <-> H+ + OH- (Kw = 1e-14)
# and acetic acid dissociation
eqsys = EqSystem.from_string("""
H2O <-> H+ + OH-; 1e-14
CH3COOH <-> CH3COO- + H+; 1.75e-5
""")

# Calculate concentrations given initial state
init_conc = {'H2O': 55.5, 'CH3COOH': 0.1, 'H+': 1e-7, 'OH-': 1e-7, 'CH3COO-': 0}
final_conc, info = eqsys.root(init_conc)

print(f"pH: {-np.log10(final_conc[eqsys.substances.index('H+')]):.2f}")

Physical Chemistry Utilities

Ionic Strength and Activity

from chempy.electrolytes import ion_strength, davies_activity_coefficient

# Calculate ionic strength of 0.1M Na2SO4
molalities = {'Na+': 0.2, 'SO4-2': 0.1}
I = ion_strength(molalities)

# Davies activity coefficient (extension of Debye-Hückel)
gamma = davies_activity_coefficient(I, z=2, eps=78.3, T=298.15)
print(f"Activity coefficient for SO4-2: {gamma:.3f}")

Practical Workflows

1. Titration Curve Simulation

def simulate_titration(acid_conc, base_concs):
    """Calculates pH as base is added to an acid."""
    results = []
    for cb in base_concs:
        # Define complex equilibrium system for each step
        eqsys = EqSystem.from_string(f"HA <-> H+ + A-; 1e-5\nH2O <-> H+ + OH-; 1e-14")
        # Solve for H+ concentration
        ...
    return results

2. Enzyme Kinetics (Michaelis-Menten ODE)

def enzyme_kinetics():
    # E + S <-> ES -> E + P
    rsys = ReactionSystem.from_string("""
    E + S <-> ES; 1e6, 1e2
    ES -> E + P; 1e3
    """)
    odesys, _ = get_odesys(rsys)
    # Integrate to see substrate depletion and product formation
    ...

3. Atmospheric Chemistry Model

def ozone_cycle():
    # Simple Chapman Cycle
    reactions = """
    O2 -> 2 O; k1
    O + O2 -> O3; k2
    O3 -> O + O2; k3
    O + O3 -> 2 O2; k4
    """
    rsys = ReactionSystem.from_string(reactions)
    # Solve for steady-state ozone concentration
    ...

Performance Optimization

Symbolic Derivation

ChemPy uses SymPy under the hood to derive the Jacobian of the ODE system. This makes integration significantly faster and more stable than numerical Jacobian estimation.

Native Code Generation

For very large systems, ChemPy can use pyodesys to generate C++ or Fortran code from your chemical reaction network, which is then compiled and called from Python.

Common Pitfalls and Solutions

The "Stiff System" Problem

Chemical systems often have reactions with widely different time scales (fast proton transfer vs. slow combustion).

# ✅ Solution: Use a stiff-capable solver (like 'cvode' or 'lsoda')
result = odesys.integrate(tout, c0, integrator='lsoda')

Formula Parsing Ambiguity

'Co' could be Cobalt or Carbon + Oxygen.

# ❌ Problem: formula_to_composition('Co') -> {27: 1} (Cobalt)
# ✅ Solution: Use Substance objects to be explicit
sub = Substance('CO', name='Carbon Monoxide')

Units and Floating Point

Equilibrium constants (K) can span 50 orders of magnitude (10⁻⁵⁰ to 10²⁰).

# ❌ Problem: Standard root finders might fail on small values
# ✅ Solution: Solve in log-space (log-concentrations)
final_conc, info = eqsys.root(init_conc, use_log=True)

ChemPy brings the precision of physical chemistry to the Python world. By automating the transition from chemical notation to mathematical models, it allows researchers to focus on the chemistry rather than the underlying differential equations.