前回に引き続きCustumTkinterを利用した軌道情報の自動計算アプリの紹介です。
今回は計算実行用の処理関数の設定コードの紹介です。
始めに入力内容から化合物の情報（charge:電荷, multi：スピン多重度）、表示軌道数（orb)を読み込んでいます。また計算する化合物群のSMILESを読み込み、split()処理を行ってリスト化を行っています。
ここで作成したSMILESのリストに従って、for文を利用して連続的に計算処理を行っています。
まずはSMILESのそれぞれの化学構造情報(smi)をもとにRdkitの機能を用いてコンフォマーサーチを行い、初期安定化構造の座標を得られるようにしています。SMILESからのデータですのでmolファイルへの変換後の初期座標はあまりあてになりませんので、このような操作を入れています。この辺りのコードについては下記の過去の記事を参照していただけると助かります。
続いて計算機条件の（threads, memory）の設定です。こちらもそれぞれの入力値をから設定できるようにしています。
計算経過の表示ですが、CustumTkinter のラベル機能でどのSMILESデータの計算が行われているか表示させるようにしています。最初は一番最初の表示のみが表示され、その後次の計算の経過は表示されなかったのですが、ウィジェットのupdate()の機能を使うことで随時の計算経過の表示ができるようになりました（ご存知の方は当たり前のことかと思いますが・・・）。
構造最適化ですが、化合物によっては所定の回数で計算が収束せずエラー（OptimizationConvergenceError）がでる場合もありますので、その対策としてtry, except での対応を行っています。（エラーで計算が途中で止まってしまうと、SMILESリストの残りの計算が行われず、未計算の部分のみを再計算を行う必要があり、非効率となるための対策です）
構造最適化の次は振動計算（処理関数：Vib_Calc())ですが、ここでは負の振動数が得られた場合（安定構造ではない）場合のチェックを行っており、負の振動数が得られた場合にはfreq_data='True' として記録し、後の結果表示で確認しやすくしています。
最後は分子軌道計算(処理関数：MO_anal())です。得られた分子軌道の情報（HOMOnp, LUMOnp)はnumpyのarrayの形で処理し、for文の処理で得られるそれぞれのSMILESに対応したデータはnumpyのvstackで蓄積するようにしています。得られたデータはsmilesのIndexを作成した上でpandasのデータフレームの形で整理しています（'FreqNega'の情報も追記）。
ここで都度データフレームを作成して、csvファイルを保存しているのは、まどろっこしいように思われるかもしれませんが、量子化学計算時に化学構造によっては先のOptimizationConvergenceError以外にもトラブルが生じて、止まってしまうことがあるので、そういった事例への対応の意味となります。他に良い方法があるのかもしれませんが、長い時間かけて実施してきた計算結果が最後の計算のエラーのために、また一からというのはなかなか厳しいものがありますので・・・
def calc_run(): global fname_smi chr=charge.get() mul=multi.get() orb=orb_input.get() with open (fname_smi) as f: readsmiles=f.read() smiles=readsmiles.split() # Frequency calculations list freq_results=[] # Smile Index list smiIndex=[] # 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() psi4.fchk(wfn, str(i)+smi+'fchk') HOMO_np=np.array(HOMOdata).reshape(1,orb+1) LUMO_np=np.array(LUMOdata).reshape(1,orb+1) if i==0: HOMO_nplist=HOMO_np LUMO_nplist=LUMO_np else: HOMO_nplist=np.vstack([HOMO_nplist, HOMO_np]) LUMO_nplist=np.vstack([LUMO_nplist, LUMO_np]) # smiles Index の作成 smiIndex.append(smi) Data_HOMO=pd.DataFrame(data=HOMO_nplist, index=smiIndex, columns=Hindex) Data_LUMO=pd.DataFrame(data=LUMO_nplist, index=smiIndex, columns=Lindex) Data_HOMO['FreqNega']=freq_results Data_LUMO['FreqNega']=freq_results Data_HOMO.to_csv('Calcd_HOMO_List.csv') Data_LUMO.to_csv('Calcd_LUNO_List.csv') Labelex_4=ctk.CTkLabel(master=root, text='Calculation was finished', font=("Times",12,"bold")) Labelex_4.place(x=50, y=410)
以上、今回は計算実行用の処理関数の紹介をさせていただきました。
次回はCustumTkinter関連のコードについて紹介したいと思います。
