Zenn
⚗️

ChemPy -とりあえず有機反応を走らせる-

に公開
Python
chemistry
kinetics
chempy
tech

イントロ

こちらの記事の動作検証で、DocumentationのExamplesに紹介されている "Chemical kinetics (system of ordinary differential equations)" のコードを、Jupyter notebookに貼り直して走らせました。
この反応式の定義部分（変数rsysを定義する部分）を差し替えて、他の反応の経時変化を描かせてみます。

題材として、こちらの文献で公開されている、酢酸エチルの加水分解反応のデータを使うこととしました。

全体としては、
AcOEt + H2O <--> AcOH + EtOH
の平衡反応ですね。
この文献のTable 3. の実験条件と、実験の結果算出した速度定数（k1 = 1.23 x 10^-4 L/mol/min, k2 = 7.09 x 10^-4 L/mol/min）を使い、Figure 2. の平衡に達するまでの経時変化曲線が再現できるか、試してみます。

進め方

お手本にしているexample codeの以下の部分、stringの内容

rsys = ReactionSystem.from_string("""2 Fe+2 + H2O2 -> 2 Fe+3 + 2 OH-; 42
2 Fe+3 + H2O2 -> 2 Fe+2 + O2 + 2 H+; 17
H+ + OH- -> H2O; 1e10
H2O -> H+ + OH-; 1e-4""")  # "[H2O]" = 1.0 (actually 55.4 at RT)

と、初期濃度を与える

c0 = defaultdict(float, {'Fe+2': 0.05, 'H2O2': 0.1, 'H2O': 1.0, 'H+': 1e-2, 'OH-': 1e-12})

の文字列を書き換えていきます。

前者は
（係数）出発物質1 + （係数）出発物質2 ... -> （係数）生成物1 + （係数）生成物2 ... ; 速度定数
という記法になっています。

化学種をいろいろな表記法で書いてみる

さて、この記法において、化学種は分子式（もしくは組成式）で書くことになっています。ここも、いろいろ試してみよう、と思い立ちました。

略号で書いてみる

有機化学者がよく書く書き方で、
AcOEt + H2O <--> AcOH + EtOHと入れてみます。

# ReactionSystem
rsys = ReactionSystem.from_string("""AcOEt + H2O -> AcOH + EtOH; 1.23e-4
AcOH + EtOH -> AcOEt + H2O; 7.09e-4""")

対応する初期濃度は

c0 = defaultdict(float, {'AcOEt': 0.5578, 'H2O': 51.63, 'EtOH': 0})

にします。
しかし、最初の反応を定義する部分のセルが、"ParseException: Expected end of text, found 'Et' (at char 3), (line:1, col:4)" で止まってしまいます。略号等の定義もしていないので、これは、まあそうでしょう。

SMILESで書いてみる

では、SMILESならどうか

# ReactionSystem
rsys = ReactionSystem.from_string("""CCOC(=O)C + [O] -> CC(=O)O + CCO; 1.23e-4
CC(=O)O + CCO -> CCOC(=O)C + [O]; 7.09e-4""")

対応する初期濃度は

c0 = defaultdict(float, {'CCOC(=O)C': 0.5578, '[O]': 51.63, 'CCO': 0})

ですね。
略号の時と同じく、最初の反応定義のセルが、"ParseException: Expected end of text, found '(' (at char 4), (line:1, col:5)"で、ダメですね。元素記号、数字以外は不可みたい。

示性式で書いてみる

出発物質、生成物がそれぞれ何かわかりやすいよう、反応定義を示性式で書いてみます。

# ReactionSystem
rsys = ReactionSystem.from_string("""CH3COOC2H5 + H2O -> CH3COOH + C2H5OH; 1.23e-4
CH3COOH + C2H5OH -> CH3COOC2H5 + H2O; 7.09e-4""")

初期濃度は

c0 = defaultdict(float, {'CH3COOC2H5': 0.5578, 'H2O': 51.63, 'C2H5OH': 0})

お、これは走りました。図1. の酢酸、酢酸エチルの経時的な濃度変化が表示されます。

図1. 化学種を示性式で書いた場合の経時変化

分子式

本来の記法である分子式で書いてみます。

# ReactionSystem
rsys = ReactionSystem.from_string("""C4H8O2 + H2O -> C2H4O2 + C2H6O; 1.23e-4
C2H4O2 + C2H6O -> C4H8O2 + H2O; 7.09e-4""")

と

c0 = defaultdict(float, {'C4H8O2': 0.5578, 'H2O': 51.63, 'C2H6O': 0})

ですね。これは想定通り動きます。当然ですが、凡例が図1と異なるだけです。

図2. 化学種を分子式で書いた場合の経時変化

わざと間違えてみる

うっかり間違えた体で、酢酸エチルからエタノールが生成するところ、プロパノールと書いちゃった、というのをやってみます。

# ReactionSystem
rsys = ReactionSystem.from_string("""C4H8O2 + H2O -> C2H4O2 + C3H8O; 1.23e-4
C2H4O2 + C3H8O -> C4H8O2 + H2O; 7.09e-4""")

"ValueError: Composition violation (1: 2) in C4H8O2 + H2O -> C2H4O2 + C3H8O"、組成が一致しないよ、と指摘されてしまいました。きっちりstoichiometry見てて偉いなぁ。

グラフの縦軸

ところで、図1.、図2. の酢酸濃度の経時曲線、600 min付近で、0.5 mol/Lくらいで平衡に達しています。ところが、原報のFigure 2、時間（横軸）は同じくらいですが、1.4 mol/Lくらいで平衡に達してますよね。再現できていないじゃん、となりそうなものですが、そもそもエステルの初期濃度0.5578 mol/Lなのに、なぜか加水分解したら酢酸根が増えてしまう？
安心してください。Figure 2. の縦軸は「酸濃度」です。つまり、加水分解して発生した酢酸と、触媒として張り込んだ塩酸を足した濃度です。原報の実験項でも、反応粗液をそのまま希釈してNaOHで滴定しています。ということで、Figure 2. の縦軸は、塩酸の初期張り大体0.9 mol/Lが乗っかった数字になっているんですね。
これを勘案すると、文献の結果をchempy.kinetics.odeで再現できている、と言えそうです。

再び略号の件

ChemPyのDocumentationのExample、Balancing stoichiometry of a chemical reactionの項の２つ目の例にあるように"balance_stoichiometry"では略号が使えるようです。
chempy.chemistry.balance_stoichiometryの引数にはsubstancesがあって、略号と組成を紐づける辞書を指定することができます。
一方で、反応オブジェクトを作るchempy.chemistryのclass Reactionでは、そういう引数がないので、今のところ、「略号は使えない」ということのように思えています。
これについて知見をお持ちの方、教えて頂ければありがたいです。

総括

今回は、chempy.ReactionSystem.from_stringを使って、文字列から反応式オブジェクトを作成し、get_odesysで連立微分方程式の反応速度式にして解く、という一連の流れをExampleのコードを元にやってみました。
また、化学種は、示性式または分子式で記述することを確認しました。大きい分子だと、おそらく分子式を採用、ということになるのでしょうが、有機化合物だとちょっと何の分子かわからなくなりそう、というのが懸念事項かな。異性化反応なんかも、区別がつかないですかね。

Discussion