pip install chempy

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

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}

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

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")

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

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())

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}")

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}")

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

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

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

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

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

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