Pythonを活用した機械学習用データ作成 その２（屈折率）

　前回までの記事で機械学習用データとして分子体積データの算出について紹介してきました。今回は分子体積データを活用した屈折率のデータ作成について紹介します。
　分子の屈折率ですが、ローレンツ・ローレンツ式を用いることで分子体積と分極率データから算出することができます（Ref　高分子論文集 , Vol. 66, No. 1, pp. 24―30 (Jan., 2009)など）。分子体積は前回までの手法で算出でき、またQM9のデータには分極率の計算データも記載されていますので、これらのデータを活用して屈折率のデータベースを作成してゆきます。

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']
qm9_pol=dfqm9['alpha']

Smiles_list=[i for i in qm9_smi]
Mol_list=[Chem.MolFromSmiles(smiles) for smiles in qm9_smi]
Polar_list=[n for n in qm9_pol]
Volume_list=[]
RefIndex_list=[]

for mol, pol in zip(Mol_list, Polar_list):
    
    mol_H = Chem.AddHs(mol)
    
    # 分子体積を関数で計算 (そのままでは実測値よりかなり大きくなるため1.5倍している)
    vdwvol=estimate_volume(mol_H)*1.5
    Volume_list.append(vdwvol)
        
    # 屈折率計算（ローレンツ　ローレンツ式から）
    vol=(1/(vdwvol*(10**-30)))
    fai=(4*3.14159/3)*vol*(pol*0.148185*10**-30)
    refindex=((1+2*fai)/(1-fai))**0.5

    RefIndex_list.append(refindex)
    
print ('calculation finished')

df_qm9_RI1=pd.DataFrame(list(zip(Smiles_list, Polar_list, Volume_list, RefIndex_list)) ,index=Smiles_list, columns=['Smiles','polarizability/a.u', 'VdwVolume/Angstrom3/mol', 'Refractive Index'])
df_qm9_RI1.to_csv('qm_RI1.csv')


print ('data_save finished')

　分子体積の計算部分は前回紹介した論文（The Journal of Organic Chemistry, 68(19), 7368--7373, 2003.）の方法を用いています。QM9からのデータの読み込みですが、今回はSmilesデータに加えて分極率（alpha)データも読み込んでいます。
　分子体積の計算値ですが、算出したデータをそのまま用いて屈折率を計算すると実測値との一致があまり良くない結果となりましたので、算出したデータを1.5倍したもので屈折率を算出しています（1.5倍程度の係数をかけるのにはそれなりの説明材料はあるのですが、説明が長くなりますのでここでは割愛させていただきます。興味がありましたらご質問などいただければと思います。）屈折率の計算は「＃屈折率計算」以下のコード部分で行っています。算出したデータをリスト化してcsvファイルで保存するのは前回までの紹介と同様です。
　実測値とは必ずしも一致しない屈折率データもあるかと思いますが、「ある程度化合物ごとの傾向がつかめる」データセットを作るという意味で参考になれば幸いです。