前回の記事で機械学習用データとして屈折率データの算出について紹介してきました。今回はさらに屈折率データを活用した誘電率のデータ作成について紹介します。
分子の誘電率ですが、オンサガーの式を用いることで屈折率、双極子モーメントデータ等から算出することができます。屈折率は前回の手法で算出でき、また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ファイルで保存するのは前回までの紹介と同様です。
実測値とは必ずしも一致しないデータも多いかと思いますが、データセットを作るという意味で参考になれば幸いです。
