Pythonを活用した機械学習用データ作成 その３（誘電率）

　前回の記事で機械学習用データとして屈折率データの算出について紹介してきました。今回はさらに屈折率データを活用した誘電率のデータ作成について紹介します。
　分子の誘電率ですが、オンサガーの式を用いることで屈折率、双極子モーメントデータ等から算出することができます。屈折率は前回の手法で算出でき、またQM9のデータには双極子モーメントの計算データも記載されていますので、これらのデータを活用して誘電率のデータベースを作成してゆきます。

import pandas as pd
import re
from  rdkit import rdBase, Chem
from rdkit.Chem import AllChem
from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem import Draw
from rdkit.Chem import Descriptors
from rdkit.ML.Descriptors import MoleculeDescriptors
from rdkit.Chem import PandasTools
import numpy as np

# QM9データ読み込み
dfqm9=pd.read_csv("gdb9_prop_smilesN1.csv")
dfqm9N=dfqm9.dropna()
qm9_smi=dfqm9N['Smiles_R']
qm9_pol=dfqm9N['alpha']
qm9_dip=dfqm9N['mu']


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]
Dip_list=[n for n in qm9_dip]
Volume_list=[]
RefIndex_list=[]
Permit_list=[]

for mol, pol, dip in zip(Mol_list, Polar_list, Dip_list):
    
    # mol = Chem.MolFromSmiles(smi)
    mol_H = Chem.AddHs(mol)
    
    # 分子体積をRdkitで計算
    try:
        AllChem.EmbedMolecule(mol_H)
        vdwvol=AllChem.ComputeMolVolume(mol_H)
    except ValueError:
        pass
        
    Volume_list.append(vdwvol)
    print (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)

    
     # 誘電率計算（オンサガー式から） 
    dip_Cm=dip*3.33*10**-30
    ref2=refindex**2
    volm=vdwvol*10**-30*6.02*10**23
    
    Equ_Uhen21=(6.02*10**23*(dip_Cm**2))/(9*(8.852*10**-12)*(1.38*10**-23)*300)
    Permittivity=(((volm+Equ_Uhen21*((ref2+2)**2))+(((volm+Equ_Uhen21*(ref2+2)**2))**2+4*2*volm*volm*ref2*ref2)**0.5))/(2*2*volm)
    Permit_list.append(Permittivity)
    print (f'Per={Permittivity}')
         
print ('calculation finished')

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

df_Sum_data.to_csv('Sum_data_Per.csv')

print ('data_save finished')

　分子体積の計算部分はRdkitで算出できる値を用いています。QM9からのデータの読み込みですが、今回はSmilesデータに加えて分極率（alpha)データ、双極子モーメント（mu)も読み込んでいます。
　分子体積の計算値ですが、今回はRdkitから算出される値をそのまま用いています（前回の屈折率の場合は係数処理を実施）。屈折率の計算は「＃屈折率計算」以下のコード部分で行っています。また誘電率の計算は「＃誘電率計算」以下のコード部分で行っています。算出したデータをリスト化してcsvファイルで保存するのは前回までの紹介と同様です。
　実測値とは必ずしも一致しないデータも多いかと思いますが、データセットを作るという意味で参考になれば幸いです。