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

化学種のプロパティ

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 の世界にもたらします。化学表記から数学的モデルへの変換を自動化することで、研究者は微分方程式の詳細ではなく、化学そのものに集中できます。