openbabelによるMOPACのインプットファイルの作成_その１_単一分子

1. openbabelによるMOPACのインプットファイルの作成

condaでopenbabelをインストールします。(google colabolatoryの人は、以下のコードでminicondaをインストーし環境設定してからconda install でopenbabelをインストールします。)

minicondaのインストール

%%bash
wget https://repo.anaconda.com/miniconda/Miniconda3-py37_4.12.0-Linux-x86_64.sh
bash Miniconda3-py37_4.12.0-Linux-x86_64.sh -b -f -p /usr/local 
conda install python=3.7 --yes

環境設定

環境設定で以下の場所を追加しておきます。

import sys
_ = (sys.path.append("/usr/local/lib/python3.7/site-packages"))

openbabelのインストール

その後、openbabelをインストールします。

conda install -c conda-forge openbabel -y

openbabelの呼び出し

openbabelの呼び出しは、以下のように行います。

import openbabel
from openbabel import pybel
# import pybel
# 昔のバージョンだと、import pybelで読み込んでいましたが、バージョンが新しくなって変わりました。環境によって変えてください。

読み込みフォーマットと出力のフォーマットの数

どれくらいの読み込みのフォーマットと出力のフォーマットがあるのかを確認します。

print(len(pybel.informats), len(pybel.outformats))

資料作成時点で、読み込みが143で出力が132となっている。

具体的なフォーマットの種類の確認

読み込みのフォーマットで具体的にどんななものを確認してみます。

pybel.informatsと入力すると、辞書型で入力形式が出力されます。

pybel.informats

pybel.outformatsと入力すると、辞書型で出力形式が出力されます。

pybel.outformats

具体的なものも確認します。全部は示さないですが、mopac形式などが読み込めることがわかります。gaussianの形式など、色々なものが読み込めて変換できます。

outputの形式も色々と出力することが可能です。mopacの形式もできますし、Chem3Dなどのフォーマットに

smilesでの単一の分子の読み込むときの具体例

単一の分子を読み込むときには、readstringを使います。
例えば、smilesで分子を読み込むときには、以下のコードでできます。

m = pybel.readstring("smiles", "C1=CC=C(C=C1)O")
m

また、最終行で読み込んだ分子を表示させると画像で出力されます。

xyz座標での分子の出力

読み込んだ分子をxyzの座標で出力してみるときには、読み込んだ分子のあとに.write(xyz)と記載して出力させます。

print(m.write("xyz"))

構造最適化を行って、xyz座標での分子の出力

構造最適化も.make3D()を使うことで行うことができます。

m.make3D(forcefield='mmff94', steps=1000)
print(m.write("xyz"))

実行すると、xyzの座標が最適化されている様子が確認できます。

MOPACの形式で出力

構造最適化した分子をmopacの形式で出力させます。

mopac_format = m.write("mop")
print(mopac_format)

一行目にPUT KEYWORDS HEREという文字、その後２行空行があって、原子と座標などが出力されているのが確認できます。

初めの3行を削除する。

最初の行は、キーワードを書く行になります。
２行目と3行目は、任意の文字のメモのために使用できる行です(空行でも問題有りません。)

今回は、次のコードを実行して最初の３行を削除して新しい３行を追加したいと思います。

mopac_format = '\n'.join(mopac_format.splitlines()[3:])

新しい3行を追加する。

新しい３行として、１行目にkeywordを加えます。
今回は2行目は、メモとしてsmilesを書き込んでおきます。
また、3行目には、メモとしてinchiを書き込んでおきます。
その後、全てを合わせてmopac_codeとして出力してみます。

mopac_keyword = "PM7 EF\n" #ここにキーワードを書く
memo1 = m.write("smiles") #任意の文字_メモ書き 。今回は分子のsmilesを書き込んだ。
memo2 = m.write("inchi") #任意の文字_メモ書き。今回は分子のinchiを書き込んだ。
mopac_code = mopac_keyword + memo1 + memo2 + mopac_format
print(mopac_code)

ファイルへの書き込み

確認できたら、phenol.datというファイルを作成し、そこに作成したmopac_codeを書き込みます。

with open("phenol.dat", mode='w') as f:
    f.write(mopac_code)

作業しているところにphenol.datファイルが作成できていれば成功になります。このファイルを使ってMOPACで計算が可能です。