前回から引き続き屈折率自動算出アプリの紹介です。
具体的なコードですが、基本的なところは以前に紹介した分子軌道情報（HOMO/LUMO)の自動算出アプリと同じです。
最初にライブラリーのインポートと各種関数の設定を行っています。
屈折率はローレンツ・ローレンツ式から分子の分極率と分子体積の情報をあれば求めることができるので、今回は分極率計算を'Pol_calc()'の関数で定義しています。分極率は基本的には分散関数を用いた計算で行うことでより正確な値が得られますが、計算コストの関係で、より簡便な計算手法（分散関数を用いないレベル）でも計算できるように設定しています。
分子体積の方はRdkitのモジュール（AllChem.ComputeMolVolume）で計算したものを用いています。ここで得られる分子体積については、ローレンツ・ローレンツ式で一般的に用いられる数値よりも小さめになっていますが、分子の構造変化に伴う「傾向」を掴むという意味では使える数値ではないかということで採用でしています。
最終的に計算された分極率と分子体積から「# 屈折率計算の部分」のコードに基づき、屈折率を計算しています。
上記説明以外の部分は分子軌道情報（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 rdkit import rdBase, Chem from rdkit.Chem import AllChem import re # 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 Pol_calc(): meth_pol=method_sv_pol.get() base=baseset_sv.get() polar = psi4.properties(meth_pol+'/'+base,properties=['polarizability'], molecule=molgeom) def calc_run(): global fname_smi chr=charge.get() mul=multi.get() with open (fname_smi) as f: readsmiles=f.read() smiles=readsmiles.split() # 振動計算の結果のリスト freq_results=[] # Smile Index のリスト smiIndex=[] # 各データのリスト polar_list=[] volume_list=[] refIndex_list=[] # 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) psi4.core.clean() # 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) # 分極率計算 Pol_calc() root.update() path = '/home/winmori/'+str(i)+smi+'-molgeom.txt' with open (path) as f: alltext= f.readlines() endgyou = len(alltext) endtxt = alltext[endgyou-1].strip() print (endtxt) polar_base=re.findall(r"[\d\.\d]+",endtxt) polar_result=float(polar_base[1]) # 分子体積計算（rdkit) try: AllChem.EmbedMolecule(mol_H) rdvdwvol=AllChem.ComputeMolVolume(mol_H) except ValueError: pass # 屈折率計算 vol=(1/(rdvdwvol*(10**-30))) fai=(4*3.14159/3)*vol*(polar_result*0.148185*10**-30) refindex=((1+2*fai)/(1-fai))**0.5 # データのまとめ refIndex_list.append(refindex) polar_list.append(polar_result) volume_list.append(rdvdwvol) smiIndex.append(smi) Data_RI=pd.DataFrame(list(zip(smiIndex, polar_list, volume_list, refIndex_list)) , columns=['Smiles','polarizability/a.u', 'VdwVolume/Angstrom3/mol', 'Refractive Index']) Data_RI['FreqNega']=freq_results Data_RI.to_csv('Calcd_RI_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のパートについて紹介したいと思います。
