分子軌道情報（HOMO/LUMO)自動算出アプリ（CustumTkinter) その3

　前回に引き続きCustumTkinterを利用した軌道情報の自動計算アプリの紹介です。
　今回は量子化学計算用の処理関数の設定コードの紹介です。構造最適化、振動計算、分子軌道の処理関数を作成しています。
　まずは構造最適化（Geom_Opt)です。CustumTkinter（コードは後の記事で紹介予定）での入力値から計算方法、基底関数などデータを読み取っています。計算手法はHF、DFT、MP２を選択できるようにしていますので、それぞれの入力値に合わせて、最適化計算できるようにif文で分岐させています。今回は自動的に多くのデータをある程度の精度で処理する必要があるので、構造最適化の最大回数（geom_maxiter）、収束しきい値（optking__g_convergence）をオプションで敢えて緩めに設定しています。
　続いて振動計算（Vib_Calc）です。構造最適化と同じく入力値から計算手法を選択して実施できるようにしています。計算結果はエネルギー（energy)と波動関数情報(wfn)に格納した上で、振動数情報をwfnから振動数情報を取り出しています（振動数は直接アウトプットされないため）。振動数のリスト（freqlist)から有用な振動数を判定で取り出していますが（freqposi, freqnega)、今回は虚数解がないかどうか(安定構造か遷移状態か）を重視して、あとで使えるよう負の振動数（freqnega)をglobal指定しています。
　最後は分子軌道（MO_anal）です。分子軌道の計算（エネルギー計算）も入力値から選択して実施できるようになっています（ただし、構造最適化、振動計算、エネルギー計算は入力値ですべて同じ手法で実施されます）。また表示させる軌道エネルギーの数も入力値（orb)に応じて変えられるようになっています。軌道エネルギーも直接表示されるわけではありませんので、波動関数情報（wfn）から適宜取り出して、HOMO、LUMOのリストとして格納しています。
　

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()
    orb=orb_input.get()
    global wfn
    global HOMOdata, LUMOdata
    global Hindex, Lindex  
    
    # Molecular orbital analysis (method in HF and MP2, function in DFT)
    if meth=='HF':        
        energy, wfn = psi4.energy(meth+'/'+base, molecule=molgeom, return_wfn=True)
        
    elif meth=='DFT':
        energy, wfn = psi4.energy(func+'/'+base, molecule=molgeom, return_wfn=True)
        
    elif meth=='MP2':
        energy, wfn = psi4.energy(meth+'/'+base, molecule=molgeom, return_wfn=True)
    
    # orbital_level : number of display level
    orbital_level=orb
    orblev=np.array(list(range(((orbital_level)+1))))

    HOMOdata=[]
    LUMOdata=[]
    Hindex=[]
    Lindex=[]
    LUMO_idx=wfn.nalpha()
    HOMO_idx=LUMO_idx-1

    for i in orblev:

        nexthomo=HOMO_idx-(i)
        nextlumo=LUMO_idx+(i)
        nhomolev=wfn.epsilon_a_subset("AO","ALL").np[nexthomo]
        nlumolev=wfn.epsilon_a_subset("AO","ALL").np[nextlumo]
    
        HOMOdata.append(nhomolev)
        LUMOdata.append(nlumolev)
        Hindex.append('homo-'+str(i))
        Lindex.append('lumo+'+str(i))
    
    return HOMOdata, LUMOdata, Hindex, Lindex   
   

　以上、今回は量子化学計算用の処理関数について紹介しました。次回は主計算の処理関数について紹介したいと思います。

分子軌道情報（HOMO/LUMO）自動算出アプリ（CustumTkinter) その２

　前回に引き続きCustumTkinterを利用した軌道情報の自動計算アプリの紹介です。
今回は処理関数の設定コードの紹介です。

# 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

# Display 'Calculation start'
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)

　まずはデータ読み込み用の関数　「data_import」を設定です。データは化学構造としてSMILES形式のものが列記されたものを想定しているので、ファイルの拡張子は'smi' または’txt'としています。指定したファイルの内容を'filename_sv’として読み込んだ上で、中身のデータを後の計算用に変数fname_smiに設定しています。
　次は計算が開始された時のスタート表示用の処理関数'start_calc_button'です。計算をスタートさせる時はCustumTkinterのButton ウィジェットで行うのですが、ボタンを押せて計算がスタートしているかどうかたまに不安になることがあるので、計算がスタートしたことを表示させる様にしてみました。tkinterのafterの機能を使うことで、計算開始のスタート表示の100msec後に実際の計算の処理関数（calc_run）がスタートするようにしています。

　今回の紹介はここまでです。次回は量子化学計算用の処理関数について紹介したいと思います。
　
　

分子軌道情報（HOMO/LUMO）自動算出アプリ（CustumTkinter) その１

　データサイエンスの活用が進むにつれて、機械学習用のデータの入手、作成の重要性が増してきているように思います。最近では学習用データとして、実験データだけでなく計算機シミュレーションの活用が進められてきています。計算データは実測データよりも比較的容易に得られますが、それでも学習用に多量のデータを準備するとなるとかなりの手間がかかってきます。このような背景から既存のSMILESデータ群のデータに対して、量子化学計算を行い、構造最適化、物性算出を行うアプリを作成してみました。

分子軌道情報算出アプリの画面

　HOMO/LUMO軌道のエネルギーは各種材料設計の際にもよく参考にされるデータかと思いますので、一つめは軌道情報の自動計算アプリを作成してみました。アプリ化に関しては、Tkinterを少しスタイリッシュにできるCustumTkinterを使ってみました。
　まずは必要なライブラリーのインストールです。

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　

アプリ化用としてcustomtkinterに加えて適宜使い分ける可能性もあるため、tkiinterもインストールしています。あとはデータ処理に必要なnumpy, pandas, 化学構造処理に必要なrdkitもインストールしています。量子化学計算はPythonコードで利用可能なこれまでのブログでも紹介してきましたpsi4を利用しています。

今回は以上になります。次回以降処理関数の設定などについて紹介してゆきたいと思います。

化学構造データ変換‐Tips その３（組成式、分子量データ）

化学構造データ変換-Tips　その３になります。
今回はsmiles形式で読み込んだデータからの化学組成式、分子量データの取り出しです。
利用データですが今回もその１の事例で使用したMolecular_Netでて提供されている水和の自由エネルギーデータ（SAMPL_MOL)のデータを利用させていただきました。
（記事の最後にコードを示しています。）

Rdkitには分子情報から機械学習等に用いる各種記述子を生成するモジュール（rdMolDescriptors）が用意されており、このモジュールを用いることで組成式や分子量の情報を取り出すことができます。
その１での事例と同様にまずは既存のcsvファイル（Mol_Net_SAMPL_Mol.csv：元ファイルからファイル名変更）読み込んだ後にsmilesが表記されている列データをリストデータ（smiles)として取り出しています。
続いて、取り出したsmilesのリストからmol形式のリスト（mols)に変換した後に、それぞれ「rdMolDescriptors.CalcMolFormula」、「rdMolDescriptors._CalcMolWt」のモジュールを用いることで組成式、分子量に変換しています。
以下はデータ追加後のデータフレームの一部ですが、組成式（chemforms)、分子量（mws)が追加されているのがわかるかと思います。

組成式、分子量追加後のデータ

以上、smilesデータからの化学組成式、分子量データの取り出しについて紹介させていただきました。

import pandas as pd
from rdkit import Chem
from rdkit.Chem import AllChem, Draw
from rdkit.Chem import rdMolDescriptors

df=pd.read_csv('Mol_Net_SAMPL_Mol.csv', encoding='shift-jis')
smiles=df['smiles']

mols=[Chem.MolFromSmiles(smi) for smi in smiles]

chemforms=[rdMolDescriptors.CalcMolFormula(mol) for mol in mols]
df['chemforms']=chemforms
mws=[rdMolDescriptors._CalcMolWt(mol) for mol in mols]
df['mws']=mws

df.to_csv('SAMPL_plus_MWetc.csv')

化学構造データ変換-Tips その２（sdfファイルの読み込み）

化学構造データ変換-Tips　その２になります。
今回はsdfファイルの読み込みです。sdfファイルは、化学構造や物性データなどを記録するためのファイル形式のひとつです。
多くの化合物の化学構造とその物性データをひとつのファイルとして保存できるため、纏まったデータを扱う際にはよく利用される形式です。
一方で、化学構造を確認しながらデータ変換などを行う際にはひとつひとつの化合物を個別の分子として読み込む必要が出てきます。
PythonライブラリーのRdkitではsdfファイルを読み込み処理するモジュール（Chem.SDMolSupplier）が用意されています。
Chem.SDMolSupplierを使うことでsdfファイルを読み込んで、分子オブジェクトのリストを作成することができます。

この記事は以下の記事等を参考に書かせていただいています。

future-chem.com

qiita.com

magattaca.hatenablog.com

RDKitにはsdfファイルを読み込む方法としてSDMolSupplierとFowardSDMolSupplierの大きく二つの方法がありますが、今回はFowardSDMolSupplierを用いて紹介しています。
ForwardSDMolSupplierモジュールの特徴は、SDファイルを一度に全て読み込むのではなく、必要に応じて一つずつ読み込むということのようです。これにより、メモリの使用量を節約できるだけでなく、大きなファイルでも高速に処理できるという利点があります。

FowardSDMolSupplierの使用例として、sdfファイルを読み込み、ひとつひとつmolオブジェクトに変換するコードを紹介したいと思います。
（関連コードは記事の最後部分）
sdfファイルですが、今回はPubChemでNaphthalene-1,3-diolで検索した結果をsdfファイルとして保存したものを用いました。
FowardSDMolSupplierでファイル名を指定してsubとして読み込んでいます。この際水素原子の情報は省かれたくないので以下のオプション説明でも出てくるようにremoveHs=Falseとしています（デフォルトなので特に記載は不要かと思いますが）。
続いて読み込んだsubをそれぞれ化合物のmol形式のリスト(mols)に変換しています。この際空データをリスト化しないようにif 文以下のコードが入っています。
さらにRdkit内でそれぞれの化合物のmolオブジェクトとして扱うためにfor 文内でMolToMolFileで処理をしています。
以上のコードでsdfファイルで提供されたファイルからRdkit で扱うmolオブジェクトへの変換を行うことができます。

因みに、FowardSDMolSupplierはSDファイルを読み込む際に、いくつかのオプションを指定することができます。これらのオプションは、分子オブジェクトの生成やプロパティの取得に影響します。

• sanitize: このオプションは、分子オブジェクトを生成する際に、化学的に正しい構造に修正するかどうかを指定します。デフォルトではTrueです。sanitizeをTrueにすると、分子の原子価や芳香性などが自動的に計算されます。sanitizeをFalseにすると、分子の構造はSDファイルに記載されたままになりますが、一部の機能が使えなくなる可能性があります。
• removeHs: このオプションは、分子オブジェクトを生成する際に、水素原子を除去するかどうかを指定します。デフォルトではFalseです。removeHsをTrueにすると、分子のサイズや計算量が減りますが、水素情報が失われる可能性があります。
• strictParsing: このオプションは、SDファイルの形式が正しいかどうかを厳密にチェックするかどうかを指定します。デフォルトではTrueです。strictParsingをTrueにすると、形式が不正なSDファイルはエラーとして扱われます。

　以上、今回はsdfファイルを読み込み、molオブジェクトに変換するコードについて紹介させていただきました。

from rdkit import Chem
from rdkit.Chem import AllChem, Draw

sup=Chem.ForwardSDMolSupplier('PubChem_compound_text_Naphthalene-1,3-diol_records.sdf',removeHs=False)

mols=[mol for mol in sup if mol is not None]

for i, mol in enumerate(mols):
    Chem.MolToMolFile(mol,f'molecule{i}.mol')
    print (i, mol) 
    
print ('File conversion was finished')
    

Psi4による量子化学計算-Tkinterアプリ14

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

github.com

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

　対応可能な範囲で機能追加などができればと思います。

Psi4による量子化学計算-Tkinterアプリ13

　引き続きPsi4の計算実施用アプリについて紹介したいと思います。

　前回までに紹介したコードを実行すると以下のような立ち上げ画面が表示されます。　

　（本シリーズの1回目にも紹介したかと思いますが）

アプリ立ち上げ画面

今回はアプリ使い方について簡単に紹介したいと思います。　

一番上の「Molecules」は計算したい分子のファイルを選択する部分で、右端の「select」を押すことで下図のようなファイル選択画面が表示されます（ファイルのあるフォルダーまでは適宜移動は必要ですが）。選択可能ファイル形式はxyz又はmolファイル形式が可能です。

ファイル選択画面

　続いて「Task」として、計算したい内容を選択します。下図では構造最適化「Geometry optimization」を選択していますが、その他振動計算や分子軌道情報、紫外可視スペクトルの計算などが選択できます。

　「Caluculation　Method」では各種計算手法、汎関数、基底関数などを選択できます。現在のコードではあまり多く選択できるようにしていませんが、必要に応じてコード内に追加することで様々な計算手法に対応はできるかと思います。

　あと「option」として計算機のスレッド数、メモリの設定と分子の電荷、多重度も設定できるようにしています。

　一番下はアウトプットファイルの名前を設定で、記入したファイル名のファイルが計算プログラムがあるディレクトリ内に作成されるようになります。

　以上の設定が終わったら、真ん中下のRunボタンを押すことで計算が開始されます。

　計算が無事終了したら、「Calculation was finished」と表示されますので、それぞれの計算タスクに応じた結果が表示されるかと思います。

　簡単なアプリですが、GitHubに公開予定ですので、良かったら使ってみてもらえればと思います。

各種設定後の画面

　今回は以上になります。次回のブログでGitHubの公開情報を紹介したいと思います。