chempy

ChemPy - 物理化学与分析化学

ChemPy 提供了系统化的方式来处理化学实体和反应。它可以自动配平复杂的氧化还原反应，求解动力学常微分方程（ODE）组，并计算平衡体系中的物种分布。

适用场景

• 配平化学方程式（包括离子反应和氧化还原反应）。
• 使用 ODE 求解器模拟化学动力学（浓度随时间的变化）。
• 计算化学平衡和物种分布（例如缓冲溶液的 pH）。
• 在化学计算中处理物理常数和单位。
• 建模离子强度和活度系数（Debye-Hückel）。
• 解析化学式并计算摩尔质量。
• 生成化学反应的 LaTeX 或 HTML 表示形式。

参考文档

官方文档：https://bjodah.github.io/chempy/
GitHub：https://github.com/bjodah/chempy
搜索模式：chempy.Reaction、chempy.ReactionSystem、chempy.balance_stoichiometry、chempy.kinetics

核心原则

化学实体

表示为字符串（'H2O'、'Fe+3'）或 Substance 对象。ChemPy 可以解析这些内容，以确定组成和电荷。

反应系统

由一组 Reaction 对象组成。此系统可以转换为数学模型（ODE），其中速率由质量作用动力学或自定义速率定律定义。

单位

ChemPy 集成了 quantities 或 pint，以确保量纲正确，这对于计算速率（例如 M⁻¹s⁻¹）至关重要。

快速参考

安装

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}

关键规则

✅ 应该做

• 使用标准记法 - 遵循标准化学记法书写化学式（'H2SO4'、'Fe+3'），以确保可靠解析。
• 定义 ReactionSystem - 对于任何复杂网络，都应将反应封装在 ReactionSystem 中，以便统一管理物质和速率。
• 检查质量／电荷平衡 - 在运行动力学模拟前，务必确认反应已配平。
• 使用单位 - 尽可能使用单位库，避免时间尺度或浓度单位出错。
• 指定速率定律 - 明确说明反应是基元反应（质量作用）还是遵循自定义速率定律。
• 向量化浓度 - 向求解器提供初始浓度时使用 NumPy 数组。

❌ 不要做

• 手动配平复杂方程式 - 交由 balance_stoichiometry 处理，尤其是氧化还原反应。
• 忽略离子强度 - 在分析化学中，请记住活度系数会随离子强度变化（使用 chempy.einstein）。
• 假设平衡快速完成 - 在动力学模型中，确保“瞬时”步骤的速率常数足够大，但又不会导致数值刚性。
• 硬编码摩尔质量 - 使用 Substance.from_formula('H2O').mass 以确保精度。

反模式（绝对不要）

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

Substance 属性

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

化学动力学（chempy.kinetics）

模拟浓度随时间的变化

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

实用工作流

1. 滴定曲线模拟

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. 酶动力学（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. 大气化学模型

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

性能优化

符号推导

ChemPy 在底层使用 SymPy 推导 ODE 系统的雅可比矩阵。这使积分速度显著快于数值估算雅可比矩阵，也更加稳定。

原生代码生成

对于非常大的系统，ChemPy 可以使用 pyodesys 将化学反应网络生成 C++ 或 Fortran 代码，随后编译并从 Python 调用。

常见问题与解决方案

“刚性系统”问题

化学系统中的反应往往具有差异很大的时间尺度（快速质子转移与缓慢燃烧）。

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

化学式解析的歧义

“Co”可能表示钴，也可能表示碳 + 氧。

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

单位与浮点数

平衡常数（K）的取值范围可能跨越 50 个数量级（10⁻⁵⁰ 到 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 将物理化学的精确性带入 Python 世界。通过自动完成从化学符号到数学模型的转换，它让研究人员能够专注于化学本身，而非底层微分方程。