Pythonを活用した機械学習用データ作成 その１（分子体積2）

　前回の記事で機械学習用のデータとしてQM9のデータを活用した分子体積データの作成について紹介させていただきました。ただ、Rdkitモジュールを利用した場合計算に少し時間がかかってしまいQM9ぐらのデータ数（１０万以上）となるとPCのスペックにもよりますが、かなり時間がかかってしまいます。今回は体積計算にRdkitモジュールを利用せずにより短時間で算出できるコードを紹介します。
　分子体積の計算は過去の論文の報告例（The Journal of Organic Chemistry, 68(19), 7368--7373, 2003）を利用しています。各原子の体積を元に結合数や環構造などを考慮しながらvan del waals 体積を算出しているものになります。

import pandas as pd
import re
from rdkit import Chem
from rdkit.Chem import AllChem

#　分子体積計算関数 (from JOC)

def estimate_volume(molecule):
    """
    Reference is below:
   　The Journal of Organic Chemistry, 68(19), 7368--7373, 2003.
    """

    #making atom volume dictionary as a look up table.
    atom_vol = {
        1:   7.24, #H
        6:  20.58, #C
        7:  15.60, #N
        8:  14.71, #O
        9:  13.31, #F
        17: 22.45, #Cl
        35: 26.52, #Br
        53: 32.52, #I
        15: 24.43, #P
        16: 24.43, #S
        33: 26.52, #As
        5:  40.48, #B
        14: 38.79, #Si
        34: 28.73, #Se
        52: 36.62  #Te
    }
    
    vol = 0

    #all atom contributions
    for atom in molecule.GetAtoms():
        vol += atom_vol[atom.GetAtomicNum()];
        
    #5.92Nb
    vol -= molecule.GetNumBonds()*5.92

    #14.7Ra, 3.8Rnr
    ring_info=molecule.GetRingInfo()
    
    rings=ring_info.NumRings()   
    atom_in_rings = ring_info.AtomRings()
    num_aromatic_ring = 0
    for ring in atom_in_rings:
        aromatic_atom_in_ring = 0
        for atom_id in ring:
            ratom = molecule.GetAtomWithIdx(atom_id)
            if ratom.GetIsAromatic():
                aromatic_atom_in_ring +=1
        if aromatic_atom_in_ring==len(ring):
            num_aromatic_ring +=1
            
    vol -=14.7*num_aromatic_ring 
    
    vol -=3.8*(rings-num_aromatic_ring)
    
    return vol;


# QM9データ読み込み
dfqm9=pd.read_csv("Mol_Net_qm9.csv")
qm9_smi=dfqm9['smiles']

Smiles_list=[i for i in qm9_smi]
Mol_list=[Chem.MolFromSmiles(smiles) for smiles in qm9_smi]

Volume_list=[]

for mol in Mol_list:
    
    mol_H = Chem.AddHs(mol)
    
    vdwvol=estimate_volume(mol_H)
        
    Volume_list.append(vdwvol)
    # print (rdvdwvol)
            
print ('calculation finished')

df_qm_Vol2=pd.DataFrame(list(zip(Smiles_list, Volume_list)),index=Smiles_list, columns=['Smiles', 'VdwVolume/Angstrom3/mol'])

df_qm_Vol2.to_csv('qm_vol2.csv')

print ('data_save finished')

　分子体積の計算は最初に関数（estimate_volume(molecule)）として定義しています。QM9データの読み込みからmolファイルのリストへの変換は前回市紹介した通りです。今回は得られたモルファイルのリストに対して分子体積計算関数を適応させて体積計算を行い、最終的に算出されたデータをリスト化した後にcsvファイルとして保存しています。
　計算時間ですが、前回のRdkitモジュールAllChem.ComputeMolVolumeでは私のPCで数時間程度かかってしましましたが、今回のコードでは数分程度で終了しました。少ないデータであればRdkitで処理した方が簡単かもしれませんが、大量のデータを処理する際の参考になれば幸いです。
　次回はこれまでの紹介記事で得られたデータをもとに、屈折率の機械学習用データを作成した例について紹介したいと思います。