【化学でPython】ChemPy：化学平衡や反応速度式をさくっと解く
はじめに
この「化学でPython」シリーズでは、化学の分野で有用な Python ライブラリを紹介しています。
今回紹介するのは、ChemPy です。
ChemPy とは？
化学反応式の係数合わせから、反応速度論のシミュレーション、平衡定数の計算まで、化学計算全般をサポートするライブラリです。
例えば、手計算では面倒な複雑な化学反応式の係数合わせや、反応過程のシミュレーションが数行のコードで書けます。
インストール
pip で簡単にインストールできます。
pip install chempy
基本的な使い方
まずは、最も基本的な使い方として、化学反応式の係数合わせ（Stoichiometry） を見てみましょう。
高校化学で習うような反応式も、Pythonを使えば一瞬で解くことができます。
ここでは、アンモニアの酸化反応（オストワルト法の途中過程）を例にします。
from chempy import balance_stoichiometry
# 反応物と生成物をセットで渡します
reactants, products = balance_stoichiometry({'NH3', 'O2'}, {'NO', 'H2O'})
print("Reactants:", reactants)
print("Products:", products)
Reactants: OrderedDict({'NH3': 4, 'O2': 5})
Products: OrderedDict({'H2O': 6, 'NO': 4})
出力結果を見ると、以下のようになっていることがわかります。
しっかりと係数が合っていますね！
実践例: 反応速度論シミュレーション
実践的な例として、「逐次反応の濃度変化シミュレーション」 をやってみます。
化学反応工学や物理化学のレポートでよく出てくる、「AがBになり、BがさらにCになる」という反応です。
この反応において、中間体 B がどのように生成・消費されるかを可視化してみましょう。
1. 反応系の定義
まず、ChemPyの ReactionSystem を使って、反応式と速度定数を定義します。
from chempy import Reaction, ReactionSystem
# 反応速度定数
k1 = 0.05 # A -> B の速度定数
k2 = 0.03 # B -> C の速度定数
# 反応の定義
# Reaction({反応物}, {生成物}, 速度定数)
r1 = Reaction({'A': 1}, {'B': 1}, k1)
r2 = Reaction({'B': 1}, {'C': 1}, k2)
# 反応系の構築
rsys = ReactionSystem([r1, r2])
print("Reaction System defined:")
for rxn in rsys.rxns:
print(rxn)
Reaction System defined:
A -> B; 0.05
B -> C; 0.03
2. シミュレーションの実行と可視化
定義した反応系に対して、初期濃度と時間を設定して計算を実行します。
ChemPyは内部で常微分方程式（ODE）を解いてくれます。
from chempy.kinetics.ode import get_odesys
import matplotlib.pyplot as plt
import numpy as np
# 初期濃度 (M = mol/L)
initial_concentrations = {'A': 1.0, 'B': 0.0, 'C': 0.0}
# 時間の設定 (0秒から100秒まで)
t_span = np.linspace(0, 100, 100)
# ODEシステム（常微分方程式系）の作成
odesys, extra = get_odesys(rsys)
# シミュレーション実行
# integrateメソッドで時間を渡すと、その時刻での濃度が返ってきます
result = odesys.integrate(t_span, initial_concentrations)
# 結果のプロット
# result.yout に濃度の時系列データが入っています
plt.figure(figsize=(8, 5))
plt.plot(t_span, result.yout[:, rsys.as_substance_index('A')], label='[A]')
plt.plot(t_span, result.yout[:, rsys.as_substance_index('B')], label='[B]')
plt.plot(t_span, result.yout[:, rsys.as_substance_index('C')], label='[C]')
plt.xlabel('Time (s)')
plt.ylabel('Concentration (M)')
plt.title('Consecutive Reaction Kinetics: A -> B -> C')
plt.legend()
plt.grid(True)
plt.show()
実行結果
グラフを見ると、以下の様子が確認できます。
- [A]の減少: 反応開始とともにAは指数関数的に減少しています。
- [B]の挙動: Aから生成されるため最初は増加しますが、Cへの分解も進むため、ある時点で極大値（ピーク）を取り、その後減少します。
- [C]の増加: 最終生成物であるCは、遅れて生成され始め、最終的には全てのAがCになります（時間が十分経てば）。
このように、ChemPyを使えば微分方程式を自分で立てて解くプログラムを書かなくても、反応式と定数を与えるだけで複雑な挙動をシミュレーションできます。
まとめ
今回は ChemPy を紹介しました。
-
Point 1: 化学反応式の係数合わせが
balance_stoichiometry一発で可能です。 - Point 2: 反応速度論のシミュレーションが、反応式を定義するだけで簡単に行えます。
- Point 3: 物理化学や化学工学の学習・研究における「検算」や「可視化」に非常に強力です。
特に、実験データの解析や、反応メカニズムの検討を行う際に、サクッとシミュレーションできるのは大きな強みです。
ぜひ試してみてください。
Discussion