Pyscfによる量子化学計算-Tkinterアプリ_11
引き続きPyscfの計算実施用アプリについて紹介したいと思います。
前回までに紹介したコードを実行すると以下のような立ち上げ画面が表示されます。
（本シリーズの1回目にも紹介したかと思いますが）
今回はアプリ使い方について簡単に紹介したいと思います。
基本的には以前に紹介したPsi4のアプリとほぼ同じ形式です。
一番上の「Molecules」は計算したい分子のファイルを選択する部分で、右端の「select」を押すことで下図のようなファイル選択画面が表示されます（ファイルのあるフォルダーまでは適宜移動は必要ですが）。選択可能ファイル形式はxyz又はmolファイル形式が可能です。
続いて「Task」として、計算したい内容を選択します。構造最適化「Geometry optimization」やその他振動計算や分子軌道情報、紫外可視スペクトルの計算などが選択できます。
「Caluculation Method」では各種計算手法、汎関数、基底関数などを選択できます。現在のコードではあまり多く選択できるようにしていませんが、必要に応じてコード内に追加することで様々な計算手法に対応はできるかと思います。
あと「option」として計算機のスレッド数、メモリの設定と分子の電荷、多重度も設定できるようにしています。
一番下はアウトプットファイルの名前を設定で、記入したファイル名のファイルが計算プログラムがあるディレクトリ内に作成されるようになります。
以上の設定が終わったら、真ん中下のRunボタンを押すことで計算が開始されます。
計算が無事終了したら、「Calculation was finished」と表示されますので、それぞれの計算タスクに応じた結果が表示されるかと思います。
簡単なアプリですが、良かったら使ってみてもらえればと思います。
Pyscfによる量子化学計算-Tkinterアプリ_10
引き続きpyscfの計算実施用アプリについて紹介したいと思います。
今回は計算実行部分の関数とTkinter設定部分のコードになります。
基本的なコードは以前に紹介したPsi4のアプリのコードと同様になりますがPyscf用に一部変更しています。計算用のインプットデータですが、Pyscfではxyzファイルでの入力に対応していますので、mol ファイルの場合はrdkitを利用して配座探索を行った後にxyz を作成し、インプットデータ(comp_geo)としています。もともとxyzファイルの場合はそのままインプットデータとしています（配座探索が現時点では組み込まれていないので留意が必要です）。以前のアプリと同様に変数によってはglobal定義にしないとエラーが出ましたのでここではglobal宣言をしていますが、他のやり方もあるのかもしれません。
def calc_run(self): global meth global func global base global svname meth=method_sv.get() func=function_sv.get() base=baseset_sv.get() svname=savename_sv.get() global chr global mul chr=charge.get() mul=multi.get() global mem global thre mem=memory.get() thre=threads.get() global fname fname=filename_sv.get() # Conformer search by rdkit then make xyz file # 計算のgeometry file(xyz)は comp_geo で設定 global comp_geo with open (fname) as f: if '.mol' in fname: mol_H=f.read() 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 = 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) with open ('geomerty_for_calc.xyz', mode='w') as f: f.write(xyz) comp_geo='geomerty_for_calc.xyz' # 上記xyxデータをファイルとして保存 # Setting input file elif '.xyz' in fname: comp_geo=fname # 計算関数内のgto(atom=###)に設定できるxyzファイルとして設定 # Select calculation task task=task_sv.get() if task=='Geometry optimization': Geom_Opt() elif task=='Vibration analysis': Vib_Calc() elif task=='Geom opt + Vib analysis': Geom_Opt_Vib_Anal() elif task=='Molecular orbitals analysis': MO_anal() elif task=='UV-Vis spectrum': UV_VIS_Spec() elif task=='Muliken charges': Mul_Charge() elif task=='Dipole moment': Dip_Moment() elif task=='Polarizability': Polar() elif task=='Infrared Spectrum': infrared_spect() # finish program def scry_finish(): exit()
以下はTkinter の設定部のコードになります。
GUIのタイトルと枠（大きさの設定）、Label, Entry, Button, Comboboxなどの設定を行っていますが、選択項目の少しの変更以外はPsi4のアプリの際と大きな変更点はありません。
Tkinter main root = tk.Tk() root.title("Pyscf Calculation Setup") root.geometry('800x650') # Select molecule file Label1=ttk.Label(root, text=u'Molecule',font=("Times","14","bold")) Label1.place(x=20, y=60) filename_sv = tk.StringVar() filenameEntry = ttk.Entry(width=60, text="", textvariable=filename_sv) filenameEntry.place(x=20, y= 90) Button1 = ttk.Button(text=u'Select',width=10) Button1.bind("<Button-1>", data_import) Button1.place(x=600, y=90) # Select calculation task Label2=ttk.Label(text=u'Task', font=("Times","14","bold")) Label2.place(x=20, y=140) task_sv = tk.StringVar() task_contents=('Geometry optimization','Vibration analysis','Geom opt + Vib analysis', 'Molecular orbitals analysis', 'UV-Vis spectrum','Muliken charges','Dipole moment', 'Polarizability', 'Infrared Spectrum' ) comboBox2=ttk.Combobox(root, height=5, width=20, state='readonly', values=task_contents, textvariable=task_sv) comboBox2.place(x=20, y=170) # Select calculation methods Label3=ttk.Label(text=u'Calculation Method', font=("Times","14","bold")) Label3.place(x=20, y=220) Label3_1=ttk.Label(root, text=u'Method',font=("Times","12")) Label3_1.place(x=20, y=240) method_sv = tk.StringVar() methods=('HF','DFT', 'MP2') comboBox3_1=ttk.Combobox(root, height=5, width=10, state='readonly', values=methods, textvariable=method_sv) comboBox3_1.place(x=20, y=260) Label3_2=ttk.Label(root, text=u'Function',font=("Times","12")) Label3_2.place(x=200, y=240) function_sv = tk.StringVar() functions=('','b3lyp','cam-b3lyp', 'edf2','m06', 'pbe','wb97x-d') comboBox3_2=ttk.Combobox(root, height=5, width=10, state='readonly', values=functions, textvariable=function_sv) comboBox3_2.place(x=200, y=260) Label3_3=ttk.Label(root, text=u'Basis set',font=("Times","12")) Label3_3.place(x=400, y=240) baseset_sv = tk.StringVar() base_sets=('3-21g','6-31g', '6-31g(d)','6-311g', 'aug-cc-pvtz') comboBox3_3=ttk.Combobox(root, height=5, width=10, state='readonly', values=base_sets, textvariable=baseset_sv) comboBox3_3.place(x=400, y=260) # Select options Label4=ttk.Label(text=u'Options', font=("Times","14","bold")) Label4.place(x=20, y=300) Label4_1=ttk.Label(root, text=u'Tread',font=("Times","12")) Label4_1.place(x=30, y=320) threads=tk.IntVar(value=2) textBox1_1=ttk.Entry(root, width=5, textvariable=threads) textBox1_1.place(x=30, y=340) Label4_2=ttk.Label(root, text=u'Memory/MB',font=("Times","12")) Label4_2.place(x=130, y=320) memory=tk.IntVar(value=500) textBox1_2=ttk.Entry(root, width=5, textvariable=memory) textBox1_2.place(x=130, y=340) Label4_3=ttk.Label(root, text=u'Charge',font=("Times","12")) Label4_3.place(x=30, y=370) charge=tk.IntVar(value=0) textBox2_1=ttk.Entry(root, width=5, textvariable=charge) textBox2_1.place(x=30, y=390) # spinのカウントは不対電子の数を数えるため通常（一重項）spin=0, ラジカルspin=1, 三重項 spin=2 Label4_4=ttk.Label(root, text=u'Multiplicity',font=("Times","12")) Label4_4.place(x=130, y=370) multi=tk.IntVar(value=0) textBox2_2=ttk.Entry(root, width=5, textvariable=multi) textBox2_2.place(x=130, y=390) # Input the name of calculated output files Label5=ttk.Label(text=u'Name of output files', font=("Times","14","bold")) Label5.place(x=20, y=470) savename_sv = tk.StringVar() textBox3=ttk.Entry(root, width=30, textvariable=savename_sv) textBox3.place(x=30, y=500) # Calculation Run Button2=ttk.Button(text=u'Run',width=20) Button2.bind("<Button-1>", calc_run) Button2.place(x=300, y=540) Label6=ttk.Label(text=u'Finish the program') Label6.place(x=630, y=560) Button3 = ttk.Button(text=u'Quit',width=10, command=scry_finish) Button3.place(x=630, y=580) root.mainloop()
以上、今回は計算実行部分とTkinterの入力部の続きについて紹介してきました。次回は実際のアプリの入力の様子の説明を行いたいと思います。
Pyscfによる量子化学計算-Tkinterアプリ_9
引き続きpyscfの計算実施用アプリについて紹介したいと思います。
今回は赤外線吸収スペクトル表示の計算の部分になります。赤外線吸収スペクトルは、これまでのブログでは紹介していませんが、Tkinterのアプリを作成するにあたり新たに追加しました。実行するにはモジュールのインポートの紹介の記事の際に記載していたかもしれませんが、以下のモジュールをインポートしておく必要があります。
from pyscf.prop.freq import rks from pyscf.prop.infrared.rks import Infrared
以下のコードはこれまでと同様に計算用の関数設定です。Tkinterアプリからの入力値（method, function, baseset等)を受けて、計算の実行を行っています。methodの違い（HF, DFT, MP2)はif文で分岐させていますが、赤外線吸収スペクトル表示については、私の方で試した範囲ではDFTのみで可能でした。そのため、HF、MP２では計算できない旨の表示をさせるようにしています。
def infrared_spect(): # Calculation of infrared spectrum # Display 'start' Label22=ttk.Label(text=u'Calculations(Infrared Spectrum) were started', font=("Times","12")) Label22.place(x=20, y=570) if meth=='HF': # HFではInfrared Spectrum 計算でエラー。計算不可の表示 Label22.place_forget() Label11_1=ttk.Label(text=u'Infrared Spectrum by HF can not be performed', font=("Times","12")) Label11_1.place(x=20, y=570) return elif meth=='DFT': mol_ir=gto.M(atom=comp_geo, basis=base,verbose=4,output=svname+'-dftir.xyz', max_memory = mem) mol_ir.spin=mul mol_ir.charge=chr lib.num_threads(thre) mf_dft_ir=dft.RKS(mol_ir) mf_dft_ir.xc=func mf_dft_ir=mf_dft_ir.newton() mf_dft_ir.chkfile=svname+'-dftir.chk' mf_dft_ir.kernel() wave_ren, modes = rks.Freq(mf_dft_ir).kernel() mf_ir_spec=Infrared(mf_dft_ir).run() elif meth=='MP2': # MP2 ではInfrared Spectrum 計算でエラー。計算不可の表示 Label22.place_forget() Label11_1=ttk.Label(text=u'Infrared Spectrum by MP2 can not be performed', font=("Times","12")) Label11_1.place(x=20, y=570) return # Display 'finish' Label22.place_forget() Label23=ttk.Label(text=u'Calculation was finished', font=("Times","12")) Label23.place(x=20, y=540)
赤外吸収スペクトルの計算は 通常のDFT計算を行った後、以下のコードを追加することで行っています。
wave_ren, modes = rks.Freq(mf_dft_ir).kernel() mf_ir_spec=Infrared(mf_dft_ir).run()
得られた計算結果は以下のコード（化学あたしいカタチさんのブログも参考）で処理することでTkinterのsubwindowに表示させています。pyscfでは計算結果に対してplot_ir()を実行するだけで簡単にスペクトル図を作成できます。
# IR Spectrum描画（化学あたらしいカタチさんの記事参考） sub_window5=tk.Toplevel() sub_window5.title('Calculations results (IR Spectrum)') sub_window5.geometry('720x540') fig=mf_ir_spec.plot_ir()[0] fig_canvas = FigureCanvasTkAgg (fig, master=sub_window5) fig_canvas.get_tk_widget().pack(fill=tk.BOTH, expand=True) toolbar = NavigationToolbar2Tk(fig_canvas, sub_window5) toolbar.update() fig_canvas.get_tk_widget().pack(fill=tk.BOTH, expand=True)
以下の図はアセトニトリルの計算で得られるスペクトル図になります。
以上、今回は赤外線吸収スペクトル表示の関数のコードについて紹介させていただきました。
次回は計算実行部分の関数とTkinter設定部分のコードについて紹介したいと思います。
Pyscfによる量子化学計算-Tkinterアプリ_8
引き続きpyscfの計算実施用アプリについて紹介したいと思います。
今回は分極率の計算の部分になります。
以下のコードはこれまでと同様に計算用の関数設定です。Tkinterアプリからの入力値（method, function, baseset等)を受けて、計算の実行を行っています。methodの違い（HF, DFT, MP2)はif文で分岐させています。分極率は、以前の記事でも紹介したようにPySCFでは、Polarizability()関数を使用します。分極率はテンソル（αxx, αxy、αxz・・・等）の値がnumpy arrayの形で算出されますので、等方的な数値のisotropic α 値を、numpyモジュールで処理をすることで算出しています。Mullken電荷や双極子モーメントの場合と同様に MP２レベルでの計算ではエラーが発生しましたので、計算できない旨のアラートを表示させるようにしています。
def Polar(): # Calculation of Polarizability # Display 'start' Label20=ttk.Label(text=u'Calculations(polarizability analysis) were started', font=("Times","12")) Label20.place(x=20, y=570) # Calculation of polarizability if meth=='HF': mol_pol=gto.M(atom=comp_geo, basis=base,verbose=4,output=svname+'-hfpol.xyz', max_memory = mem) mol_pol.spin=mul mol_pol.charge=chr lib.num_threads(thre) mf_hf_pol=scf.RHF(mol_pol) mf_hf_pol.chkfile=svname+'-hfpol.chk' mf_hf_pol.kernel() polar=mf_hf_pol.Polarizability().polarizability() elif meth=='DFT': mol_pol=gto.M(atom=comp_geo, basis=base,verbose=4,output=svname+'-dftpol.xyz', max_memory = mem) mol_pol.spin=mul mol_pol.charge=chr lib.num_threads(thre) mf_dft_pol=dft.RKS(mol_pol) mf_dft_pol.xc=func mf_dft_pol=mf_dft_pol.newton() mf_dft_pol.chkfile=svname+'-dftpol.chk' mf_dft_pol.kernel() polar=mf_dft_pol.Polarizability().polarizability() elif meth=='MP2': # MP2 ではPolarizability 計算でエラー。計算不可の表示 Label20.place_forget() Label11_1=ttk.Label(text=u'Polarizability calculations by MP2 can not be performed', font=("Times","12")) Label11_1.place(x=20, y=570) return # Display 'finish' Label20.place_forget() Label21=ttk.Label(text=u'Calculation was finished', font=("Times","12")) Label21.place(x=20, y=540)
以上、今回は分極率の計算の関数のコードについて紹介させていただきました。
次回は赤外吸収スペクトル表示の関数のコードについて紹介したいと思います。
Pyscfによる量子化学計算-Tkinterアプリ_7
引き続きpyscfの計算実施用アプリについて紹介したいと思います。
今回は双極子モーメントの計算の部分になります。
以下のコードはこれまでと同様に計算用の関数設定です。Tkinterアプリからの入力値（method, function, baseset等)を受けて、計算の実行を行っています。methodの違い（HF, DFT, MP2)はif文で分岐させています。双極子モーメントは、以前の記事でも紹介したようにdip_moment()の関数を用いることで計算可能です。計算結果（dipole)はベクトル量であり、各成分（x, y, z方向）とその大きさが出力（Debye単位）されます。 また分子全体の双極モーメントの大きさ(dipole_magnitude )はxyzの各成分の2乗和のルートで求められるので、numpyで以下の処理を行うことで算出できます。Mullken電荷の場合と同様に MP２レベルでの計算ではエラーが発生しましたので、計算できない旨のアラートを表示させるようにしています。
def Dip_Moment(): #Calculation of Dopole moment # Display 'start' Label17=ttk.Label(text=u'Calculations(Dipole analysis) were started', font=("Times","12")) Label17.place(x=20, y=570) # Molecular orbital analysis (method in HF and MP2, function in DFT) if meth=='HF': mol_dip=gto.M(atom=comp_geo, basis=base,verbose=4,output=svname+'-hfdip.xyz', max_memory = mem) mol_dip.spin=mul mol_dip.charge=chr lib.num_threads(thre) mf_hf_dip=scf.RHF(mol_dip) mf_hf_dip.chkfile=svname+'-hfdip.chk' mf_hf_dip.kernel() dipole_vec=mf_hf_dip.dip_moment() elif meth=='DFT': mol_dip=gto.M(atom=comp_geo, basis=base,verbose=4,output=svname+'-dftdip.xyz', max_memory = mem) mol_dip.spin=mul mol_dip.charge=chr lib.num_threads(thre) mf_dft_dip=dft.RKS(mol_dip) mf_dft_dip.xc=func mf_dft_dip=mf_dft_dip.newton() mf_dft_dip.chkfile=svname+'-dftdip.chk' mf_dft_dip.kernel() dipole_vec=mf_dft_dip.dip_moment() elif meth=='MP2': # MP2 ではdipole moment 計算でエラー。計算不可の表示 Label17.place_forget() Label11_1=ttk.Label(text=u'Dipole moment calculations by MP2 can not be performed', font=("Times","12")) Label11_1.place(x=20, y=570) return # Display 'finish' Label17.place_forget() Label18=ttk.Label(text=u'Calculation was finished', font=("Times","12")) Label18.place(x=20, y=540) # Calculation of dipole moment # Display results of dipole monment from x, y, x vector data Dipmom=round(np.sqrt(np.sum(dipole_vec **2)),3) Label19=ttk.Label(text=f'Dipole Moment= {Dipmom} D', font=("Times","12")) Label19.place(x=20, y=570)
以上、今回は双極子モーメントの電荷計算の関数のコードについて紹介させていただきました。
次回は分極率計算の関数のコードについて紹介したいと思います。
Pyscfによる量子化学計算-Tkinterアプリ_6
今回はマリケンの電荷計算の部分になります。
以下のコードはこれまでと同様に計算用の関数設定です。Tkinterアプリからの入力値（method, function, baseset等)を受けて、計算の実行を行っています。methodの違い（HF, DFT, MP2)はif文で分岐させています。Muliken 電荷は、以前の記事でも紹介したようにmulliken_pop()の関数を用いることで計算可能です。計算結果はmul_popに格納されます。MP２レベルでの計算（そもそもMullken電荷はあまり高レベルでの計算には適してはいませんが）ではエラーが発生しましたので、計算できない旨のアラートを表示させるようにしています。
def Mul_Charge(): # Calculation of Mulliken Charges # Display 'start' Label15=ttk.Label(text=u'Calculations(Mulliken analysis) were started', font=("Times","12")) Label15.place(x=20, y=570) # Molecular orbital analysis (method in HF and MP2, function in DFT) if meth=='HF': mol_mul=gto.M(atom=comp_geo, basis=base,verbose=4,output=svname+'-hfmul.xyz', max_memory = mem) mol_mul.spin=mul mol_mul.charge=chr lib.num_threads(thre) mf_hf_mul=scf.RHF(mol_mul) mf_hf_mul.chkfile=svname+'-hfmul.chk' mf_hf_mul.kernel() mul_pop=mf_hf_mul.mulliken_pop(mol_mul) elif meth=='DFT': mol_mul=gto.M(atom=comp_geo, basis=base,verbose=4,output=svname+'-dftmul.xyz', max_memory = mem) mol_mul.spin=mul mol_mul.charge=chr lib.num_threads(thre) mf_dft_mul=dft.RKS(mol_mul) mf_dft_mul.xc=func mf_dft_mul=mf_dft_mul.newton() mf_dft_mul.chkfile=svname+'-dftmul.chk' mf_dft_mul.kernel() mul_pop=mf_dft_mul.mulliken_pop(mol_mul) elif meth=='MP2': # MP2 ではMulliken Charge 計算でエラー。計算不可の表示 Label15.place_forget() Label11_1=ttk.Label(text=u'Mulliken Charge calculations by MP2 can not be performed', font=("Times","12")) Label11_1.place(x=20, y=570) return # Display 'finish' Label15.place_forget() Label16=ttk.Label(text=u'Calculation was finished', font=("Times","12")) Label16.place(x=20, y=570)
以下は電荷データを表示させるコードです。電荷に関するデータは第二要素に格納されているため、そのデータを「Charge」として取り出しています。またそれらデータを各原子に対応させたるためfor文以下の処理を行うと共に、pandasのデータフレームとして表示させる形にしています。表示についてはこれまでと同様にTkinterのToplevelのサブウィンドウを活用する形です。また、以前にも紹介させていただいた化学あたらしいカタチさんの記事を参考に化学構造上に電荷をマッピングした図を作成し、別の画面で表示させるようにもしています（sub_window4=tk.Toplevel()以下）。
# Calculation of Mulliken charges charge=mul_pop[1] atom_symbol=[] mulliken_data=[] for n in range(len(mol_mul._atom)): print (mol_mul.atom_symbol(n), round(charge[n],4)) atom_symbol.append(mol_mul.atom_symbol(n)) mulliken_data.append(round(charge[n],4)) Mullikendf=pd.DataFrame(({'Atom': atom_symbol, 'Mulliken Charge':mulliken_data})) # Display results of Mullken data sub_window3=tk.Toplevel() sub_window3.title('Calculations results (Mulliken Data)') sub_window3.geometry('720x540') LabelS_9=ttk.Label(sub_window3, text='Mulliken Charges', font=("Arial","14", 'bold')) LabelS_9.place(x=10, y=10, width=300) LabelS_10=ttk.Label(sub_window3, text=Mullikendf, font=("Arial","12")) LabelS_10.place(x=10, y=50) # 電荷の画像描画（化学あたらしいカタチさんの記事参考） sub_window4=tk.Toplevel() sub_window4.title('Calculations results (Mulliken Charges)') sub_window4.geometry('720x540') mol_file1=Chem.MolFromXYZFile(comp_geo) rdDetermineBonds.DetermineBonds(mol_file1) print (Chem.MolToMolBlock(mol_file1)) mol_file1H=Chem.AddHs(mol_file1) fig=SimilarityMaps.GetSimilarityMapFromWeights(mol_file1H, mul_pop[1], colorMap='RdBu') fig.savefig('fig1.png', bbox_inches='tight') canvas1 = tk.Canvas(sub_window4, bg="white", height=500, width=700) canvas1.pack() figure1=tk.PhotoImage(file='fig1.png', height=500, width=700) canvas1.create_image(0, 0, image=figure1, anchor=tk.NW) sub_window4.mainloop()
以上、今回はマリケンの電荷計算の関数のコードについて紹介させていただきました。
次回は双極子モーメント計算の関数のコードについて紹介したいと思います。
Pyscfによる量子化学計算-Tkinterアプリ_5
引き続きpyscfの計算実施用アプリについて紹介したいと思います。
今回は紫外可視吸収スペクトルの表示の部分になります。スペクトル表示に関しては以前の記事でも紹介した下記のHPのコードを参考にしました。ここでは吸収強度は算出された遷移エネルギーを中心に一定のGauss分布を示すとして計算しています。
https://github.com/jamesETsmith/2022_simons_collab_pyscf_workshop/blob/main/demos/05_Excited_States.ipynb
# setting of spectrum analysis def gaussian (x, mu, sig): return np.exp(-np.power(x-mu, 2.)/(2*np.power(sig, 2.))) def spectral_analysis(mf_tddft): spectrum_width=0.1 osc_strengths=mf_tddft.oscillator_strength()[:states] print (osc_strengths) energies_ev=mf_tddft.e[:states]*ha_2_ev x_range=np.linspace(energies_ev.min()*0.9, energies_ev.max()*1.1, num=1000) intensity=np.zeros(x_range.size) for e, f in zip(energies_ev, osc_strengths): intensity +=gaussian(x_range, e, spectrum_width)*f dx=(x_range[-1]-x_range[0])/x_range.size area=(intensity*dx).sum() intensity /=area return x_range, intensity,osc_strengths, energies_ev
以下のコードはこれまでと同様に計算用の関数設定です。Tkinterアプリからの入力値（method, function, baseset等)を受けて、計算の実行を行っています。methodの違い（HF, DFT, MP2)はif文で分岐させています。今回は紫外可視吸収スペクトルの表示ということでTD計算（tddft.TDDFT())の実行を行っています。methodの違い（HF, DFT, MP2)はif文で分岐させています。計算させるstate数（励起状態に相当）はここでは「10」を設定しています。
def UV_VIS_Spec(): global states states=10 # Display 'start' Label14=ttk.Label(text=u'Calculations(UV_VIS_Spectrum) were started', font=("Times","12")) Label14.place(x=20, y=570) # TDSCF calculation (method in HF and MP2, function in DFT) if meth=='HF': mol_u=gto.M(atom=comp_geo, basis=base,verbose=4,output=svname+'-hfuv.xyz', max_memory = mem) mol_u.spin=mul mol_u.charge=chr lib.num_threads(thre) mf_hf_u=scf.RHF(mol_u) mf_hf_u.chkfile=svname+'-hfuv.chk' mf_hf_u.kernel() mf_tddft=tddft.TDDFT(mf_hf_u) mf_tddft.nstates=states mf_tddft.kernel() mf_tddft.analyze() elif meth=='DFT': mol_u=gto.M(atom=comp_geo, basis=base,verbose=4,output=svname+'-dftuv.xyz', max_memory = mem) mol_u.spin=mul mol_u.charge=chr lib.num_threads(thre) mf_dft_u=dft.RKS(mol_u) mf_dft_u.xc=func mf_dft_u=mf_dft_u.newton() mf_dft_u.chkfile=svname+'-dftuv.chk' mf_dft_u.kernel() mf_tddft=tddft.TDDFT(mf_dft_u) mf_tddft.nstates=states mf_tddft.kernel() mf_tddft.analyze() elif meth=='MP2': mol_u=gto.M(atom=comp_geo, basis=base,verbose=4,output=svname+'-mp2uv.xyz', max_memory = mem) mol_u.spin=mul mol_u.charge=chr lib.num_threads(thre) mf_mp_uN=mol_u.RHF().run() mf_mp2_u=mf_mp_uN.MP2() mf_mp2_u.chkfile=svname+'-mp2uv.chk' mf_mp2_u.kernel() mf_tddft=tddft.TDDFT(mf_mp2_u) mf_tddft.nstates=states mf_tddft.kernel() mf_tddft.analyze() # Display 'finish' Label14.place_forget() Label15=ttk.Label(text=u'Calculation was finished', font=("Times","12")) Label15.place(x=20, y=570)
下記からはスペクトルの表示画面の設定です。以前のPsi4の場合と同様にTkinterのtoplevelのsub_windowを活用して表示させています。Sub_window1はデータフレームにまとめた波長、振動子強度の結果の表示です。スペクトル表示は別の画面（sub_window2)で行っています。計算データから得られた結果をmatplotlibでグラフ作成していますが、Tkinter上で表示させるためにFigurecanvasTkAggを活用しています。またグラフ表示調整を画面上で行えるようにするためNavigationToolbar2Tkも活用しています。 この辺も以前のPsi4のTkinterアプリと同様です。
uvdata = {"Excitation Energy (eV)":[], "Intensity":[], "Excitation Energy (nm)":[]}
x_range, intensity, osc_strengths, energies_ev=spectral_analysis(mf_tddft)
waverength=1240/x_range
uvdata["Excitation Energy (eV)"] +=x_range.tolist()
uvdata["Intensity"] +=intensity.tolist()
uvdata["Excitation Energy (nm)"] += waverength.tolist()
energies_nm0=np.array([1240/energies_ev])
energies_nm0=np.ravel(energies_nm0)
print (energies_nm0)
print (osc_strengths)
osc_strengths_N=np.array([0 if e=="NaN" else e for e in osc_strengths])
# osc_strengths_N.reshape(1,10)
print (osc_strengths_N)
print (energies_nm0.shape)
print (osc_strengths_N.shape)
Ex_nm_strength={'Waverength/nm':energies_nm0, 'Os_Strength':osc_strengths_N}
print (Ex_nm_strength)
nm_strength_df=pd.DataFrame(Ex_nm_strength)
# Display energy data and spectrum
sub_window1=tk.Toplevel()
sub_window1.title('Calculations results (Energy Data)')
sub_window1.geometry('720x540')
LabelS_7=ttk.Label(sub_window1, text='Excitation Energy and Oscilleator Strength', font=("Arial","14", 'bold'))
LabelS_7.place(x=10, y=10, width=300)
LabelS_8=ttk.Label(sub_window1, text=nm_strength_df, font=("Arial","12"))
LabelS_8.place(x=10, y=50)
sub_window2=tk.Toplevel()
sub_window2.title('Calculations results (Spectrum)')
sub_window2.geometry('720x540')
fig = Figure()
ax = fig.add_subplot(1,1,1)
fig, ax = pltls.subplots()
x_data=uvdata['Excitation Energy (nm)']
y_data=uvdata['Intensity']
ax.plot(x_data, y_data)
pltls.xlabel('$\lambda$ / nm')
pltls.ylabel('Intensity / unit')
fig_canvas = FigureCanvasTkAgg (fig, master=sub_window2)
fig_canvas.get_tk_widget().pack(fill=tk.BOTH, expand=True)
toolbar = NavigationToolbar2Tk(fig_canvas, sub_window2)
toolbar.update()
fig_canvas.get_tk_widget().pack(fill=tk.BOTH, expand=True)
pltls.show()
以上、今回は分紫外可視吸収スペクトルの表示の部分コードについて紹介させていただきました。
次回はマリケンの電荷計算の関数 のコードを紹介してゆきたいと思います。
