chempy

# ChemPy - Physical and Analytical Chemistry

ChemPy provides a systematic way to handle chemical entities and reactions. It can automatically balance complex redox reactions, solve systems of ordinary differential equations (ODEs) for kinetics, and calculate species distribution in equilibria.

## When to Use

- Balancing chemical equations (including ionic and redox reactions).
- Simulating chemical kinetics (concentration vs. time) using ODE solvers.
- Calculating chemical equilibria and speciation (e.g., pH of a buffer).
- Handling physical constants and units in chemical calculations.
- Modeling ionic strength and activity coefficients (Debye-Hückel).
- Parsing chemical formulas and calculating molar masses.
- Generating LaTeX or HTML representations of chemical reactions.

## Reference Documentation

**Official docs**: https://bjodah.github.io/chempy/  
**GitHub**: https://github.com/bjodah/chempy  
**Search patterns**: `chempy.Reaction`, `chempy.ReactionSystem`, `chempy.balance_stoichiometry`, `chempy.kinetics`

## 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

```bash
pip install chempy
```

### Standard Imports

```python
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

```python
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)

```python
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)

```python
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

```python
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

```python
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

```python
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

```python
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

```python
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)

```python
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

```python
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).

```python
# ✅ 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.

```python
# ❌ 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²⁰).

```python
# ❌ 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.