RDKitにおける、原子と結合の扱い方

## はじめに
**RDKitを使った基本的な原子と結合の取り扱い**について解説します。

※マテリアルズインフォマティクス関係の内容を他にも投稿していますので、よろしければ[こちら](https://qiita.com/oki_kosuke/items/34b1071224953776c211#%E6%8A%95%E7%A8%BF%E4%B8%80%E8%A6%A7)の一覧から他の投稿も見て頂けますと幸いです。

### 環境

- windows10
- conda 4.10.3
- python 3.7.10
- rdkit 2021.03.5

## 原子の取り扱い
例として、トルエンを描きます。

```python
from rdkit import Chem
# toluene
Chem.MolFromSmiles('Cc1ccccc1')
```

![toluene.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1057114/ed8ca2f7-1181-8224-18ca-52fe274aca6f.png)

化合物を図として保存したい場合は以下のコードで保存します。

```python
# 画像の保存
from rdkit.Chem import Draw
toluene = Chem.MolFromSmiles('Cc1ccccc1')
Draw.MolToFile(toluene, 'toluene.png')
```

### 水素原子の取り扱い
通常、RDKitでは水素原子は暗に扱われ表には現れません。例として、トルエン（C<sub>7</sub>H<sub>8</sub>）の原子数を数えてみます。

```python
toluene.GetNumAtoms()
```
```:出力
7
```

ベンゼン環と側鎖のメチル基の炭素原子の合計は7であり、それに付随している水素原子は原子数にカウントされていないことが分かります。立体構造を考える際など水素原子を表したい場合は以下のコードで水素原子を付加することができます。

```python
toluene_h2 = Chem.AddHs(toluene)
# 原子数を確認
toluene_h2.GetNumAtoms()
```
```:出力
15
```

先ほどはカウントされていなかった水素原子もカウントされていることが分かります。水素原子が付加されたことは以下のように構造式を描いても確認することができます。

```python
toluene_h2
```

![toluene_h2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1057114/2cfc9b0b-9bce-6335-700e-8408c288bf45.png)

また、一度追加した水素原子は以下のコードで除去することができます。

```python
# 水素原子を除去
toluene_remove_h2 = Chem.RemoveHs(toluene_h2)
# 原子数をカウント
toluene_remove_h2.GetNumAtoms()
```
```:出力
7
```

### 化合物中の各原子の取り扱い
例として以下の化合物を考えます。適当な構造を描いてSMILESを取得しました。

```python
m = Chem.MolFromSmiles("NC(C1=CC(C(CN2)C3C2=O)=C3C=N1)/C=C/C4CO4")
m
```

![m.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1057114/d293d404-2fe9-c63f-7046-523f01c4a922.png)

以下の議論でどの原子に着目しているか分かるように、各原子のインデックスを表示します。

```python
# 各原子に番号を振っておく
m_with_num = Chem.MolFromSmiles("NC(C1=CC(C(CN2)C3C2=O)=C3C=N1)/C=C/C4CO4")
for atom in m_with_num.GetAtoms():
    atom.SetProp("atomLabel", str(atom.GetIdx()))
m_with_num
```

![m_with_num.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1057114/cdc20acc-9083-79d5-f6ff-eb51c9c4c8b4.png)

注目している原子の原子番号は以下のコードで確認できます。

```python
# 原子番号の取り出し（0番目の原子の原子番号を表示）
m.GetAtoms()[0].GetAtomicNum()
```
```:出力
7
```

0番目の原子の原子番号は7（＝窒素N）ということが分かります。

### 原子と環の関係性
着目している原子が環に含まれるかは以下のコードで確認できます。

```python
# 各原子が環に含まれるか判定する（例：6番の原子が環に含まれるか判定する）
m.GetAtomWithIdx(6).IsInRing()
```
```:出力
True
```

6番目の原子は環に含まれているという判定です。実際に6番目の原子の炭素は右側の5員環に含まれています。

以下のコードで化合物中の各原子が環に含まれているかを全て判定できます。

```python
# 原子数を定義しておく
num = m.GetNumAtoms()
# 繰り返し処理で「対象原子の番号、環に含まれるか否か」を出力する
for i in range(num):
    print(str(i)+"番目の原子", m.GetAtomWithIdx(i).IsInRing())
```
```:出力
0番目の原子 False
1番目の原子 False
2番目の原子 True
3番目の原子 True
4番目の原子 True
5番目の原子 True
6番目の原子 True
7番目の原子 True
8番目の原子 True
9番目の原子 True
10番目の原子 False
11番目の原子 True
12番目の原子 True
13番目の原子 True
14番目の原子 False
15番目の原子 False
16番目の原子 True
17番目の原子 True
18番目の原子 True
```

再度、検討している化合物を示します（インデックスつき）

![m_with_num.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1057114/f9b45527-6e69-cc43-6936-1cff53d9bdbb.png)

また、以下のコードで着目している原子がある数の環に含まれるか判定できます。

```python
# 4番目の原子が4員環に含まれるか判定する
print(m.GetAtomWithIdx(4).IsInRingSize(4))
# 4番目の原子が5員環に含まれるか判定する
print(m.GetAtomWithIdx(4).IsInRingSize(5))
# 4番目の原子が6員環に含まれるか判定する
print(m.GetAtomWithIdx(4).IsInRingSize(6))
```
```:出力
True
False
True
```

次のコードで化合物中のどの原子がどの数の環に含まれるのか表示することができます。

```python
# どの環にどの番号の原子が含まれるのか表示する
for i in (3,4,5,6):
    print(str(i)+"員環に含まれる原子の番号")
    for j in range(num):
        if m.GetAtomWithIdx(j).IsInRingSize(i) == True:
            print(j)
```
```:出力
3員環に含まれる原子の番号
16
17
18
4員環に含まれる原子の番号
4
5
8
11
5員環に含まれる原子の番号
5
6
7
8
9
6員環に含まれる原子の番号
2
3
4
11
12
13
```

## 結合の取り扱い
続いて、結合の取り扱いについて解説致します。まずは例で示した化合物の結合にインデックスを付けます。[こちら](https://rdkit.blogspot.com/2020/04/new-drawing-options-in-202003-release.html)を参考にしました。

```:python
from rdkit.Chem.Draw import rdMolDraw2D
from IPython.display import SVG

d2d = rdMolDraw2D.MolDraw2DSVG(600,300)
d2d.drawOptions().addBondIndices=True
d2d.DrawMolecule(m)
d2d.FinishDrawing()
SVG(d2d.GetDrawingText())
```

![m_with_bondnum.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1057114/57611037-c052-0d44-d772-da5af997601a.png)

### 結合の種類
各結合の結合様式を表示することができます。

```python
# 0番目の原子の結合タイプを表示
print(m.GetBonds()[0].GetBondType())
# 2番目の原子の結合タイプを表示
print(m.GetBonds()[2].GetBondType())
# 9番目の原子の結合タイプを表示
print(m.GetBonds()[9].GetBondType())
```
```:出力
SINGLE
AROMATIC
DOUBLE
```

各結合のインデックス、その結合の開始と終了の原子のインデックス、結合のタイプを表示します。

```python
# 各結合を形成する原子の番号とその結合の対応を表示
for bond in m.GetBonds():
    print(str(bond.GetIdx())+"番目の結合", "開始原子："+str(bond.GetBeginAtomIdx()), "終了原子："+str(bond.GetEndAtomIdx()), "　結合のタイプ：", bond.GetBondType())
```
```:出力
0番目の結合 開始原子：0 終了原子：1 　結合のタイプ： SINGLE
1番目の結合 開始原子：1 終了原子：2 　結合のタイプ： SINGLE
2番目の結合 開始原子：2 終了原子：3 　結合のタイプ： AROMATIC
3番目の結合 開始原子：3 終了原子：4 　結合のタイプ： AROMATIC
4番目の結合 開始原子：4 終了原子：5 　結合のタイプ： SINGLE
5番目の結合 開始原子：5 終了原子：6 　結合のタイプ： SINGLE
6番目の結合 開始原子：6 終了原子：7 　結合のタイプ： SINGLE
7番目の結合 開始原子：5 終了原子：8 　結合のタイプ： SINGLE
8番目の結合 開始原子：8 終了原子：9 　結合のタイプ： SINGLE
9番目の結合 開始原子：9 終了原子：10 　結合のタイプ： DOUBLE
10番目の結合 開始原子：4 終了原子：11 　結合のタイプ： AROMATIC
11番目の結合 開始原子：11 終了原子：12 　結合のタイプ： AROMATIC
12番目の結合 開始原子：12 終了原子：13 　結合のタイプ： AROMATIC
13番目の結合 開始原子：1 終了原子：14 　結合のタイプ： SINGLE
14番目の結合 開始原子：14 終了原子：15 　結合のタイプ： DOUBLE
15番目の結合 開始原子：15 終了原子：16 　結合のタイプ： SINGLE
16番目の結合 開始原子：16 終了原子：17 　結合のタイプ： SINGLE
17番目の結合 開始原子：17 終了原子：18 　結合のタイプ： SINGLE
18番目の結合 開始原子：13 終了原子：2 　結合のタイプ： AROMATIC
19番目の結合 開始原子：9 終了原子：7 　結合のタイプ： SINGLE
20番目の結合 開始原子：11 終了原子：8 　結合のタイプ： SINGLE
21番目の結合 開始原子：18 終了原子：16 　結合のタイプ： SINGLE
```

## まとめ
RDKitを使った基本的な原子と結合の取り扱いについて解説いたしました。

### 参考
- [RDKit blog](https://rdkit.blogspot.com/)
- [RDKit日本語ドキュメント](https://www.rdkit.org/docs_jp/Getting_Started_with_RDKit_in_Python_jp.html)