Zenn
🔥

RDKitで作った分子から各元素の数を取り出して燃やしてみる

2022/06/23に公開
Python
RDKit
tech

はじめに

分子式から、炭素の数、酸素の数、といった情報を抜き出したいと思いました。
例えば、完全燃焼の反応式を記述したり、その際に必要な最小の酸素量＝理論酸素量を算出するためには、こういった情報が必要になります。

C_xH_yO_z + \frac{(2x + 0.5y - z)}{2} O_2　-> xCO_2 + \frac{y}{2}H_2O　\tag{式1}

式1であれば、分子中の各原子の数x, y, zから理論酸素量 (2x + 0.5y - z) / 2 を求めることができます。
ところが意外なことに、rdkitで、入力した分子のx, y, zを、個別に求める機能を見つけることができませんでした。
よく使われているH社の化学描画ソフトCと表計算ソフトM社のEを組み合わせると、それ用の関数CHEM_NUM_ATOMSが準備してあり、手軽に処理できます。
rdkitは、pythonで化学構造を扱う際に非常に有用なのですが、同じことを単独では行う方法が見つからず、他のライブラリと組み合わせることにより実現できました。

概要

rdkitで分子を構成する全原子の種類をリストとして取得し、標準ライブラリであるcollections.Counter()で種類毎の出現回数を数えることで、分子中の各原子の数を数えることができた。

内容

1. 事例用の分子を作成し、分子式を取得したものの…

化合物として、（気分で決めた）プロピオン酸エチルのmolオブジェクトmを作ります。分子式は、ChemモジュールのサブモジュールrdMolDescriptorsの関数を使い、以下で取得できます。

from rdkit import Chem
from rdkit.Chem import rdMolDescriptors 
m = Chem.MolFromSmiles('CCC(=O)OCC')
rdMolDescriptors.CalcMolFormula(m)

[出力]: 'C5H10O2'

これは文字列ですので、元素記号と、それに並ぶ数字を別々に取得できれば、目的を達成できます。ところが、元素記号は一文字のこともあれば二文字のこともあります。ある元素の数も、一桁のこともあれば、二桁以上のこともあります。このため、文字数で切り分ける、といった手段も使えず、良いロジックを思いつけなかったため、この方法は断念しました。

2. 課題解決の方針を設定した

rdkitでは、分子を構成する各原子の種々の属性を取得することができます。この属性の中には、「原子の種類＝元素」も含まれています。今回は、この属性を利用し、以下のロジックを用いることにより、課題を解決することができました。

  1. 分子を構成する全ての原子の々の種類をリストとして取得する。
  2. 要素の種類（今回の場合は元素）毎の出現回数を数える。
  3. 要素をkey、出現回数をvalueとする辞書型のリストにまとめる。
  4. keyを指定して、valueを呼び出すことにより、分子式中のそれぞれの原子の数を特定する。

3. 水素原子を含む各原子を取得した

以下のコードで、項目1で作成したプロピオン酸エチルのmolオブジェクトmを構成する原子の種類を個々に取得し、リストにしました。

#失敗例  
l =[]
for A in m.GetAtoms():
    a=A.GetSymbol()
    l.append(a)
print(l)

[出力]: ['C', 'C', 'C', 'O', 'O', 'C', 'C']

しかし、このリストには、「水素原子」が含まれていません。水素原子の情報も含めて、分子を構成する原子をリストの要素とするためには、原子の種類の取得に先立ち、元となるmolオブジェクトに、下記の操作Chem.AddHsで水素を付加しておく必要があります。

m2 = Chem.AddHs(m)
l =[]
for A in m2.GetAtoms():
    a=A.GetSymbol()
    l.append(a)
print(l)

[出力]: ['C', 'C', 'C', 'O', 'O', 'C', 'C', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H']

これで、分子を構成する原子毎に、原子の種類を取得し、リストとすることができました。

4. リストの要素の種類毎に出現回数を数えた

pythonの標準ライブラリであるcollectionsモジュールの、サブクラスCounterを使うと、リスト中の要素の出現回数を、要素の種類をkey、出現回数をvalueとする辞書型のリストとして取り出すことができます。

import collections
d = {}
d = collections.Counter(l) 
print(d)
d.get('C', 0)

Counter({'H': 10, 'C': 5, 'O': 2})
[出力]: 5

次いで、getメソッドで、元素'C'の出現回数を取り出しています。なお、辞書型のリストから、[ ]内にkeyを指定してvalueを取り出すこともできますが、存在しないkeyを指定してしまうと、エラーが発生します。このため、その元素が存在しない場合には、出現回数として0を返すように指定して、getメソッドを使いました（補足参照）。今回の例では、「炭素の数」を求めているので、分子式にもあるように、「5」が正しく返されています。

5. 理論酸素量を求めてみた

最後に、式1の理論酸素量を求めてみました。適用できる分子を拡げるため、構成元素としてC, H, Oの他、N, S, Pも加えました（有機化合物のかなりの割合はカバーできだろう、と）。なお、N, S, Pは完全燃焼して、それぞれNO_2, SO_2, P_2O_5となる、としました。
各元素の数に、その元素が完全燃焼するのに必要な酸素原子数をかけ、計算用のリストtaoに追記した後、リストの要素の総和を求め、酸素分子数とするために、2で割りました。また、酸素との量論関係を指定していない元素が分子式に含まれる場合には、燃焼の反応式を定義するよう、警告を出すようにしています（tryからexceptのブロック）。

TAO = 0
tao = []
for k in d.keys():
    nO = 0
    v = d[k]
    if k == 'C' or k == 'N' or k == 'S':
        nO = v*2
    elif k == 'H':
        nO = v*0.5
    elif k == 'O':
        nO = -v
    elif k == 'P':
        nO = v*2.5
    else:
        nO = 'xx'
    tao.append(nO)
try:
    TAO = sum(tao)/2
    print('Theoretical Amount of Oxygen (O2) = ', TAO)
except:
    print('define reaction formula')

[出力]: Theoretical Amount of Oxygen (O2) = 6.5

プロピオン酸エチルの分子式中には、N, S, Pはありませんが、正しい理論酸素量6.5（酸素分子として）が返されています。もちろん、N, S, Pを含む分子を入力しても、正しい理論酸素量が返されるおことは、別途確認しました。　　
分子式から、完全燃焼に最小限必要な酸素量を求める、という、化学者にとってはなんということもない計算ですが、コードにすると30行程度を要しました。

まとめ

rdkitに該当する機能がない場合でも、pythonの他のライブラリや関数と組み合わせることにより、目的を達成できました。
pythonは汎用的に活用できるプログラミング言語で、rdkitのように化学に特化したライブラリや、scikit-learnのような、回帰、分類といったデータ解析のためのライブラリ以外にも、色々な、有用な機能を実現できるライブラリがあります。自身がやりたいことに応じて、種々のライブラリの使い方を学んでいくのも大事だと思いました。また、同じことを、複数の方法で実現できる場合もあるため、「どういう方法、ロジックが最適か」を想像＋創造していく能力も重要だと思います。

補足

python3.8では、

d['P']

と分子中に存在しない元素をkeyとして指定しても、valueとして0が返ってきましたので、本文の記法と同じ結果となり、どちらの書き方でも良いことになります。ただし、公式チュートリアルには「存在しないキーを使って値を取り出そうとするとエラーになります。」とあります。また、python3.7以前で、この記法で存在しないkeyを指定すると、エラーが返ってくるため、本文中の記法としました。

Discussion