前回から引き続き紫外可視吸収波長自動算出アプリの紹介です。
具体的なコードですが、基本的なところは以前に紹介した分子軌道情報（HOMO/LUMO)や屈折率の自動算出アプリと同様です。
最初にライブラリーのインポートと各種関数の設定を行っています。
紫外可視吸収波長（エネルギー）と振動子強度の計算には「tdscf_excitations」をインポートしています。励起エネルギーの計算や波長データへの変換等は以前のブログ記事「Psi4による量子化学計算-紫外可視光吸収データ（その１，その２）」を参照してもらえればと思います。
今回のアプリでは設定項目として、計算する励起エネルギーレベル数として「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のパートについて紹介したいと思います。
