chempy

What this skill does

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