Pythonで立体配座（コンフォーマー）を生成する

## はじめに
分子は結合軸周りの回転によってさまざまな３次元構造をとる。Pythonで分子の３次元構造を生成する方法を紹介する。ここで生成した３次元構造はGaussianなどの量子化学計算の良い初期構造になりうる。

## 手順
①３次元構造を生成
②分子力学法によって最適化
今回のプログラムではこれらを一連の動作として実行

### 環境
- python  3.7.13
- rdkit  2022.3.4
- py3dmol  1.8.1
- pandas  1.2.5

### 参考サイト
- [３次元構造の生成](https://future-chem.com/rdkit-3dstruct)
- [RDkitを用いた分子の操作（3D構造の最適化）](https://qiita.com/ottantachinque/items/80a22d2428231bcdfc74#%E4%B8%89%E6%AC%A1%E5%85%83%E6%A7%8B%E9%80%A0%E3%81%AE%E6%9C%80%E9%81%A9%E5%8C%96%E3%81%AB%E3%81%A4%E3%81%84%E3%81%A6)
- [分子力学法による最適化](https://future-chem.com/rdkit-conformer)

## 3次元構造の生成

3次元構造生成には
- ディスタンス・ジオメトリー（DG）法
- KDG
- ETDG
- ETKDG\
およびその改良法があるが、**基本的にはETKDGを使えばよい**。

DG法に対して、
- ベンゼン環は平面
- アルキンは直線 \
といった一般的な化学常識を加えたのがKDGである。またケンブリッジ結晶構造データベースの実験データから得た二面角の分布傾向（experimental torsion angle preferences）を考慮したものがETDGであり、それらの両方を加えたのがETKDGである。

```Python
#ライブラリーのインポート
from rdkit import rdBase, Chem
from rdkit.Chem import AllChem
from rdkit.Chem.Draw import IPythonConsole
import py3Dmol
import pandas as pd
```

データの読み込み。ここでは例としてSMILESが分かっている4つの分子が記載されたcsvファイルを読み込みます。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1057114/f1fe6aaa-e837-09b9-5ba9-1efe7845b67c.png)

```Python
# データの取得
df = pd.read_csv("sample_products.csv")
df
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1057114/cad03773-4579-7289-bf2f-289a8a1bf889.png)

各化合物のSMILESをmol形式に変換しておきます。また水素原子を付加した状態のmol形式の取得しておきます。RDKitの公式にも記述がありますが「良い３Dコンフォメーションを得るには、ほとんどの場合、まず初めに水素原子を分子に付加すること上手くいきます。

```python
# データフレームの化合物名を定義（のちのファイル生成で使う）
labels = df["ID"]
# 各化合物のmol（水素付加なし）を取得
mols = [Chem.MolFromSmiles(mol) for mol in df["SMILES"]]
# 各化合物のmol（水素付加済み）を取得
molHs = [ Chem.AddHs(Chem.MolFromSmiles(mol)) for mol in df["SMILES"]]
```

コンフォーマー生成のパラメータを定義しておきます。AllChem.EmbedMultipleConfsの引数として代表的なものは以下です。

- どの分子の3次元構造を生成するか（mol）
- いくつのコンフォマーを生成するか（numConfs）
- どの程度類似した構造を排除するか（pruneRmsThresh）

RMSの値を大きくすると得られるコンフォーマーの数は減少します。


```python
# コンフォーマー生成のパラメータの定義
confs = 1000
rms = 0.5
```

## フォースフィールド（力場）をかけて最適化する

以下のコードでは

1. 1000個のコンフォマーを発生
2. 各コンフォマーをUFFまたはMMFFで構造最適化（UFFもしくはMMFFの力場をかける）
3. 各コンフォマーのエネルギーを計算
4. 最安定構造から順番にソートして（相対エネルギー，コンフォマーID）のタプルのリストを取得

という作業を行っています。

```python
# コンフォーマーの生成とフォースフィールド法で構造最適化
def opt_sp_mm(m_h, confs, rms, ff):
    # 1. 1000個のコンフォマーを発生（ETKDG法）
    cids = AllChem.EmbedMultipleConfs(m_h, numConfs=confs, randomSeed=1234, pruneRmsThresh=rms, numThreads=0)

    # 2,3. 各コンフォマーを最適化し，エネルギー計算
    energy = []
    if ff == 'uff':
        for cid in cids:
            uff = AllChem.UFFGetMoleculeForceField(m_h, confId=cid)
            uff.Minimize()
            energy.append((uff.CalcEnergy(), cid))
    if ff == 'mmff':
        prop = AllChem.MMFFGetMoleculeProperties(m_h)
        for cid in cids:
            mmff = AllChem.MMFFGetMoleculeForceField(m_h, prop, confId=cid)
            mmff.Minimize()
            energy.append((mmff.CalcEnergy(), cid))

    # 4. エネルギーをソートし，相対エネルギーとIDをリストに格納
    energy.sort()
    return [(i-energy[0][0],j) for i,j in energy]
```
まずはuffで最適化してエネルギー計算を実行し、最安定構造から順番に5つのコンフォーメーションを抽出します。この後Gaussianで量子化学計算を行うことを想定して、Gaussianの計算実行ファイルを作るためのテキストを準備します。

```python
# uffとmmffで各化合物の最適構造を求めていく
uff_gaussian_filenames = []

# 最適化した構造の上位何件を取得するか設定する
pickup = 5

for i in range(len(labels)):    
    # uff法での構造最適化
    uff_e = opt_sp_mm(molHs[i], confs, rms, 'uff')
    # uffで最適化した構造のID上位を取得する
    uff_confIds = [m for l,m in uff_e[:pickup]]
    
    # 繰り返し回数の定義
    num_confIds = len(uff_confIds)
    if num_confIds<pickup:
        iter = num_confIds
    else:
        iter = pickup
    # 上位をmolファイルで保存していく
    for j in range(iter):
        uff_filename = f"{labels[i]}_uff_"+str(j)+".mol"
        uff_gaussian_filename = f"g16 {labels[i]}_uff_"+str(j)+".gjf"
        uff_gaussian_filenames.append(uff_gaussian_filename)
        Chem.MolToMolFile(molHs[i], confId=uff_confIds[j], filename=uff_filename)    
```

Gaussianの計算実行ファイルを作るためのテキストを作成します。

```python
# Gaussian計算実行時のテキストを作成
txt = ["#!/bin/bash"]
txt.extend(uff_gaussian_filenames)
# 各要素に改行コードを付加
obj = map(lambda x: x + "\n", txt)
with open("uff_gaussian.txt", "w", encoding="utf-8") as f:
    f.writelines(obj)
```
次にmmffでも同様の手順でコンフォーメーションを抽出していきます。
```python
# uffとmmffで各化合物の最適構造を求めていく
mmff_gaussian_filenames = []

# 最適化した構造の上位何件を取得するか設定する（DC-02がconfs = 1000、rms = 0.5で3件しかヒットしない）
pickup = 5

for i in range(len(labels)):    
    # mmff法での構造最適化
    mmff_e = opt_sp_mm(molHs[i], confs, rms, 'mmff')
    # muffで最適化した構造のID上位を取得する
    mmff_confIds = [m for l,m in mmff_e[:pickup]]
    
    # 繰り返し回数の定義
    num_confIds = len(mmff_confIds)
    if num_confIds<pickup:
        iter = num_confIds
    else:
        iter = pickup
    # 上位をmolファイルで保存していく
    for j in range(iter):
        mmff_filename = f"{labels[i]}_mmff_"+str(j)+".mol"
        mmff_gaussian_filename = f"g16 {labels[i]}_mmff_"+str(j)+".gjf"
        mmff_gaussian_filenames.append(mmff_gaussian_filename)
        Chem.MolToMolFile(molHs[i], confId=mmff_confIds[j], filename=mmff_filename)
```
```python
# Gaussian計算実行時のテキストを作成
txt = ["#!/bin/bash"]
txt.extend(mmff_gaussian_filenames)
# 各要素に改行コードを付加
obj = map(lambda x: x + "\n", txt)
with open("mmff_gaussian.txt", "w", encoding="utf-8") as f:
    f.writelines(obj)
```

個人的な感覚としてはuffとmmffであまり大きな差異はない印象。ただ、力場をかけた方が何もしない（ETKDG法で求めた）構造よりも量子化学計算後の最適化エネルギーがより安定な気がする。