あるケミストの独り言(winchemwinの日記)

ケミスト(化学者)の視点で、面白そうな情報(シミュレーション関係など)を発信

紫外可視吸収波長自動算出アプリ(CustumTkinter)その2

前回から引き続き紫外可視吸収波長自動算出アプリの紹介です。
 具体的なコードですが、基本的なところは以前に紹介した分子軌道情報(HOMO/LUMO)や屈折率の自動算出アプリと同様です。
 最初にライブラリーのインポートと各種関数の設定を行っています。
 紫外可視吸収波長(エネルギー)と振動子強度の計算には「tdscf_excitations」をインポートしています。励起エネルギーの計算や波長データへの変換等は以前のブログ記事「Psi4による量子化学計算-紫外可視光吸収データ(その1,その2)」を参照してもらえればと思います。
 今回のアプリでは設定項目として、計算する励起エネルギーレベル数として「states」を設定しています(デフォルトは10)。
 Stete数に応じて得られたデータはnumpyのアレイ型データとした後に、SMILES、振動計算の結果とともにpandas のデータフレーム型としてまとめています。各列に同じSMILES、振動計算の結果が並ぶお世辞にも見栄えの良いデータとはなっていませんが、データの確認には支障がないかと思ってます。また時間があればこの辺も改良できればと思います。
 
 上記説明以外の部分は分子軌道情報(HOMO/LUMO)の自動算出アプリとほぼ同じかと思います。

import customtkinter as ctk
import tkinter as tk
from tkinter import ttk
from tkinter import StringVar, messagebox
from turtle import onclick
from tkinter import filedialog
import numpy as np
import pandas as pd
from numpy import matlib
from tkinter.ttk import Entry
import psi4
from psi4.driver.procrouting.response.scf_response import tdscf_excitations
from  rdkit import rdBase, Chem
from rdkit.Chem import AllChem


# Smiles file select by filedialog 
def data_import(self):
    type=[('smi file','*.smi'), ('text file', '*.txt')]
    file=filedialog.askopenfilename(filetypes=type)
    filename_sv.set(file)
    global fname_smi
    fname_smi=filename_sv.get()

    return fname_smi

def start_calc_button(self):
    Labelex_1=ctk.CTkLabel(master=root, text='Calculation was started',font=("Times",12))
    Labelex_1.place(x=50, y=370)
    
    root.after(100, calc_run)
    
def Geom_Opt():
    # Geometry optimization
    # getting meth(calculation method), func(function), base(basis set) 
    meth=method_sv.get()
    func=function_sv.get()
    base=baseset_sv.get()
    
    psi4.set_options({'geom_maxiter': 50,
    'optking__g_convergence': 'gau_loose',})
    
    # geometry optimize calculation (method in HF and MP2, function in DFT)
    if meth=='HF':
        psi4.optimize(meth+'/'+base, molecule=molgeom)
        
    elif meth=='DFT':
        
        psi4.optimize(func+'/'+base, molecule=molgeom)
        
    elif meth=='MP2':
        psi4.optimize(meth+'/'+base, molecule=molgeom)
    
    optimized_geom = molgeom.save_string_xyz_file()
                  
def Vib_Calc():
    # Frequency calculation
    meth=method_sv.get()
    func=function_sv.get()
    base=baseset_sv.get()
    
    # Frequency calculation (method in HF and MP2, function in DFT)
    if meth=='HF':        
        energy, wfn = psi4.frequency(meth+'/'+base, molecule=molgeom, return_wfn=True)
        
    elif meth=='DFT':
        energy, wfn = psi4.frequency(func+'/'+base, molecule=molgeom, return_wfn=True)
        
    elif meth=='MP2':
        energy, wfn = psi4.frequency(meth+'/'+base, molecule=molgeom, return_wfn=True)
    
    # making frequency list

    freqlist = np.array(wfn.frequencies())
    freqposi=freqlist[(freqlist>1)&(freqlist<5000)]
    global freqnega
    freqnega=freqlist[(freqlist<0)]

    return freqposi, freqnega

def MO_anal():
    # Molecular orbital analysis
    meth=method_sv.get()
    func=function_sv.get()
    base=baseset_sv.get()
    sta=states_input.get()
    global wfn
    global ene_mo 
    global sumarray
    
    psi4.set_options({'save_jk':True})
    
    # Molecular orbital analysis (method in HF and MP2, function in DFT)
    if meth=='HF':        
        ene_mo, wfn = psi4.energy(meth+'/'+base, molecule=molgeom, return_wfn=True)
    
    elif meth=='DFT':
        ene_mo, wfn = psi4.energy(func+'/'+base, molecule=molgeom, return_wfn=True)
        
    elif meth=='MP2':
        ene_mo, wfn = psi4.energy(meth+'/'+base, molecule=molgeom, return_wfn=True)
    
    # TDDFT calculation (UV-VIS)  
    res = tdscf_excitations (wfn, states =sta)

    exenergy = [r["EXCITATION ENERGY"] for r in res]
    osstrength = [r["OSCILLATOR STRENGTH (LEN)"] for r in res]

    arrayex = np.array(exenergy)
    arrayos = np.array(osstrength)
    arraynm0 =[(1240/(arrayex*27.21162))]
    arraynm = np.array(arraynm0)

    arrayexr = np.reshape(arrayex,(-1,1))
    arrayosr = np.reshape(arrayos,(-1,1))
    arraynmr = np.reshape(arraynm,(-1,1))
    
    sumarray=np.c_[arrayexr,arraynmr,arrayosr]
    

def calc_run():
    global fname_smi
    chr=charge.get()
    mul=multi.get()
    sta=states_input.get()
    
    with open (fname_smi) as f:
        readsmiles=f.read()
    
    smiles=readsmiles.split()
    
    # Frequency calculations list
    freq_results=[]
    
    #  Conformer search by rdkit then make xyz file
    for i, smi in enumerate(smiles): 

        mol = Chem.MolFromSmiles(smi)
        mol_H = Chem.AddHs(mol)
            
        confs = AllChem.EmbedMultipleConfs(mol_H, 10, pruneRmsThresh=1)
        prop = AllChem.MMFFGetMoleculeProperties(mol_H)

        energy=[]
        
        for conf in confs:
            mmff = AllChem.MMFFGetMoleculeForceField(mol_H, prop,confId=conf)
            mmff.Minimize()
            energy.append((mmff.CalcEnergy(), conf))
    
            conflist = np.array(energy)
            sortconf=conflist[np.argsort(conflist[:,0]) ]

            stconfid=sortconf[0,1]
            stconfid2=int(stconfid)
            stconfgeom=mol_H.GetConformer(stconfid2)

        xyz = f'{chr} {mul}'
        for atom, (x,y,z) in zip(mol_H.GetAtoms(), stconfgeom.GetPositions()):
            xyz += '\n'
            xyz += '{}\t{}\t{}\t{}'.format(atom.GetSymbol(), x, y, z)
                  
        #  Setting input file
        global molgeom
        molgeom = psi4.geometry(xyz)
               
    # Set calculation method
        thre=threads.get()
        mem=memory.get()
        psi4.set_num_threads(nthread=thre)
        psi4.set_memory(mem*1000000)
    
        # Set output files
        psi4.set_output_file(str(i)+smi+'-molgeom.txt')

        # 計算経過の表示
        # global root, Labelex_2
        Labelex_2=ctk.CTkLabel(master=root, text='No'+str(i)+'  '+smi+' Calc was started',font=("Times",10))
        Labelex_2.place(x=50, y=390)

        Labelex_2.update()
        
        # 構造最適化ー収束エラー時の対応も含む
        try:
            Geom_Opt()
            root.update()
        except  psi4.driver.p4util.exceptions.OptimizationConvergenceError as err:
            print (err)
       
            continue
        
        optimized_geom = molgeom.save_string_xyz_file()
        with open (str(i)+smi+'-optgeom.xyz','w') as f:
            f.write(optimized_geom)
        
        Vib_Calc()
        root.update()
        if freqnega:
            freq_data='True'         
        else:
            freq_data='False'
            
        # freq_results.append(freq_data)  
          
        MO_anal()
        root.update()
        
        smiIndex=np.repeat(smi, sta).reshape(-1,1) 
        freq_result=np.repeat(freq_data, sta).reshape(-1,1)
        
        UVdata_np=np.hstack([smiIndex, sumarray])
        UVdata_np=np.hstack([UVdata_np, freq_result])
        
        if i==0:
            sumUV_nplist=UVdata_np

        else:
            sumUV_nplist=np.vstack([sumUV_nplist, UVdata_np])
            
        
        # column の作成
        columns1=["SMILES","Excitation Energy /au ","Excitation Energy /nm","Oscilleator Strength ","Freq_result"]
     
        Data_UV_VIS=pd.DataFrame(data=sumUV_nplist, columns=columns1)          
            
        Data_UV_VIS.to_csv('Calcd_UV_VIS_List.csv')     
    
    Labelex_4=ctk.CTkLabel(master=root, text='Calculation was finished', font=("Times",12,"bold"))
    Labelex_4.place(x=50, y=410)
                
# finish program]

def scry_finish():
    exit() 

 以上、紫外可視吸収波長自動算出アプリのコード(関数設定、吸収波長計算等)について紹介させていただきました。
 次回はCustumTkinterのパートについて紹介したいと思います。

紫外可視吸収波長自動算出アプリ(CustumTkinter)その1


 以前の記事でHOMO-LUMO軌道、屈折率の自動算出アプリを紹介してきました。今回以降別の物性値として紫外可視吸収波長の自動算出アプリを紹介してゆきたいと思います。基本的なコードはこれまで紹介してきたアプリと同様ですが、紫外可視吸収データを算出させるためのTDDFT計算をおこなっているところが少し異なります。皆さんの参考になれば幸いです。
 アプリの起動画面は以下のような感じになります。今回もCustumTkinterで作成しています。

紫外可視吸収算出アプリの立ち上げ画面


このアプリでは化学物質のSMILESリストから紫外可視吸収波長及び振動子強度のリストを自動作成し、csvファイルに保存できるようになっています。画面的には前述のように、これまでのアプリとほぼ同じですが、紫外可視吸収波長計算用の計算の設定の部分が少し変わっています。
 次回以降具体的なコードについて紹介してゆきたいと思います。

Pythonを活用した機械学習用データ作成 その3(誘電率)

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

Pythonを活用した機械学習用データ作成 その2(屈折率)

 前回までの記事で機械学習用データとして分子体積データの算出について紹介してきました。今回は分子体積データを活用した屈折率のデータ作成について紹介します。
 分子の屈折率ですが、ローレンツローレンツ式を用いることで分子体積と分極率データから算出することができます(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ファイルで保存するのは前回までの紹介と同様です。
 実測値とは必ずしも一致しない屈折率データもあるかと思いますが、「ある程度化合物ごとの傾向がつかめる」データセットを作るという意味で参考になれば幸いです。

Pythonを活用した機械学習用データ作成 その1(分子体積2)

 前回の記事で機械学習用のデータとしてQM9のデータを活用した分子体積データの作成について紹介させていただきました。ただ、Rdkitモジュールを利用した場合計算に少し時間がかかってしまいQM9ぐらのデータ数(10万以上)となると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で処理した方が簡単かもしれませんが、大量のデータを処理する際の参考になれば幸いです。
 次回はこれまでの紹介記事で得られたデータをもとに、屈折率の機械学習用データを作成した例について紹介したいと思います。

Pythonを活用した機械学習用データ作成 その1(分子体積1)

 機械学習を活用した材料開発が活発になってきています。機械学習には元になるデータセットが必要ですが、実験データ数には限界もあることから、計算機シュミレーションを活用したデータセットも作成されています。その中でもQMデータセットは数多くの機械学習研究においてベンチマークとして用いられている代表的なデータセットになります。QMデータセットの中でもQM9は10万以上の低分子有機化合物の高レベル量子化学計算のデータが収集されているデータセットで様々な場面で活用されています。ただQM9に収集されていないパラメーターを機械学習に使いたい場面もしばしば遭遇します。そこでQM9のデータを活用して新たに機械学習用のパラメータデータを作成するコードを作成してみましたので紹介したいと思います。今回は分子体積になります(追加できそうなパラメーターが新たにできましたら、随時紹介したいと思います。)

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

# 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)
    
    # 分子体積をRdkitで計算
    try:
        AllChem.EmbedMolecule(mol_H)
        rdvdwvol=AllChem.ComputeMolVolume(mol_H)
    except ValueError:
        pass
        
    Volume_list.append(rdvdwvol)
    # print (rdvdwvol)
            
print ('calculation finished')

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

df_qm_Vol1.to_csv('qm_vol1.csv')

print ('data_save finished')

 最初にQM9のデータセット(MoleculeNet等から入手可能)をpandasのデータフレーム(dfqm9)として読み込んでいます。分子体積の計算はRdkitライブラリーの関数でmolファイルから算出可能ですので、QM9のデータセットからmolファイルデータを作成します。
 まず化学構造情報のSMILESデータを取り出した後(qm9_smi)、molファイルに変換しています。分子体積は水素原子を加えた形から計算されますので、molファイルに水素原子を加えた(Chem.AddHs(mol))後にAllChem.EmbedMoleculeで3D構造を発生させています。得られた構造対してAllChem.ComputeMolVolumeの処理を行うことで分子体積を発生させています。最終的に算出されたデータをリスト化したものをデータフレームとした後にcsvファイルとして保存しています。
 以上の処理で分子体積は算出されますが、データ処理にある程度(私のPCでは数時間程度)かかってしまいました。上記の方法でも問題はないかと思いますが、数分程度で終了する別の方法もありますので、次回はその方法を紹介したいと思います。

屈折率自動算出アプリ(CustumTkinter) その5

 前回までに屈折率の自動計算のアプリのコード、利用画面の説明などについて紹介してきました。作成したコードをGitHubに公開していますので、MITライセンスのもとご自由に利用いただければと思います。
 

github.com

 ご利用にあたって、コメント、感想などがあればご連絡ください。