openbabel / Drkcore

openbabel-2.3.1が出てますね。

Sat, 28 Jan 2012 09:10:47 +0900

2.3.1がリリースされたようです。

個人的に興味があるのは

PNG files from Open Babel contain molecular information and can be read to give the MDL Molfile.
Pybel now uses the built-in 2D depiction, and no longer needs OASA.

とABINITのフォーマットに対応したあたりかな。

あと、openbabel-python.iをいじってたので、ここをいじった場合のコンパイルのオプションをメモっておく。swigが有効になるようにしないといけないのに気付かなくてハマった。

cmake ../openbabel-2.3.1 -DPYTHON_BINDINGS=ON -DEIGEN2_INCLUDE_DIR=/usr/local/tmp/eigen-eigen-2.0.12 -DRUN_SWIG=ON

OBGenericからOBOrbitalDataへのキャストをできるようにして、vectorの設定もしたので、手元のpythonバインディングでは

orb = toOrbitalData(mol.GetData(openbabel.ElectronicData))
orb.GetAlphaOrbitals()[orb.GetAlphaHOMO()-1].GetEnergy()

とやるとHOMOのエネルギー(eV)を得られるようになっている。

追記12.01.28

homebrewでいれたpythonで使いたい場合optionで指示する

$ cmake ../openbabel-2.3.1 -DPYTHON_BINDINGS=ON \
  -DPYTHON_LIBRARY=/usr/local/lib/libpython2.7.dylib -DPYTHON_EXECUTABLE=/usr/local/bin/python \
  -DEIGEN2_INCLUDE_DIR=/Users/kzfm/openbabel/eigen-eigen-2.0.17 -DRUN_SWIG=ON

Python用のGAMESSラッパーを書いている

Wed, 22 Jun 2011 19:07:02 +0900

去年書いたGAMESSラッパーに手を加えてGitHubにあげた。ヘッダーの生成まわりはもっとやらないといけないんだけど、基底関数とコントロールまわりは動くようにした。といっても一点計算と最適化ぐらいしかしないんだけど。

エラー終了しているときにはエラーの内容を出力できるようにした
rungmsのパスの確認
gamess inputを出力できるようにした
SCF計算のタイプも指定できる

こんな感じで動かします。例としてEthane。デフォルトはCPUに優しいSTO3Gの一点計算です。

import gamess
g = gamess.Gamess()
obc = ob.OBConversion()
obc.SetInFormat("mol")

mol = ob.OBMol()
next = obc.ReadFile(mol, "examples/ethane.mol")
print g.gamess_input(mol)
try:
    newmol = g.run(mol)
except GamessError, gerr:
    print gerr.value

print newmol.GetEnergy()
print [(obatom.GetIdx(), obatom.GetType(), obatom.GetPartialCharge()) for \
obatom in ob.OBMolAtomIter(newmol)]

結果はこれ。

 $contrl runtyp=energy scftyp=rhf  $end
 $basis gbasis=sto ngauss=3 $end
 $SYSTEM MWORDS=30 $END
 $DATA
6324
C1
C      6.0     -0.7560000000    0.0000000000    0.0000000000 
C      6.0      0.7560000000    0.0000000000    0.0000000000 
H      1.0     -1.1404000000    0.6586000000    0.7845000000 
H      1.0     -1.1404000000    0.3501000000   -0.9626000000 
H      1.0     -1.1405000000   -1.0087000000    0.1781000000 
H      1.0      1.1404000000   -0.3501000000    0.9626000000 
H      1.0      1.1405000000    1.0087000000   -0.1781000000 
H      1.0      1.1404000000   -0.6586000000   -0.7845000000 
 $END

-78.30530748
[(1, 'C3', -0.16967199999999999), (2, 'C3', -0.16967199999999999), \
(3, 'HC', 0.056557999999999997), (4, 'HC', 0.056559999999999999), \
(5, 'HC', 0.056554), (6, 'HC', 0.056559999999999999), \
(7, 'HC', 0.056554), (8, 'HC', 0.056557999999999997)]

ラジカルの計算がしたいのでROHFかUHFの設定ができるようにしておいたがスピン多重度の指定が出来ないのでとっととやる。

あと、テスト書かなあかんなぁと思いながらエキスパートPythonを読んでます。二周目か三周目かわからんけど、何回読んでもこの本は楽しい。

エキスパートPythonプログラミング
Tarek Ziade
アスキー・メディアワークス / 3780円 ( 2010-05-28 )

水素を引きぬいてラジカルにする

Sat, 18 Jun 2011 15:08:48 +0900

共有結合先の原子をチェックして冗長なものが出ないようにした。GAMESSのテストも終わったのであとは走らせるだけ。

共有結合先の原子を知るメソッドがなくて、総当りでGetBondを行う。もし結合が形成されてなければNoneが返ってくるのでそれで評価できるらしい。

import openbabel as ob
import sys

def clone(mol):
    obc = ob.OBConversion()
    obc.SetInAndOutFormats("mol", "mol")
    molstring = obc.WriteString(mol)
    new_mol = ob.OBMol()
    obc.ReadString(new_mol, molstring)
    return new_mol

def abs_h(mol):
    radicals = []
    mol.AddHydrogens()
    h_indexs = []
    hetero_indexs = []

    for atom in ob.OBMolAtomIter(mol):
        if atom.GetType() == 'H':
            h_indexs.append(atom.GetIdx())
        else:
            hetero_indexs.append(atom.GetIdx())

    hatoms = []
    non_redundant_hs = []
    for h_index in h_indexs:
        hetero = neighbor(mol, h_index, hetero_indexs)
        if hetero not in hatoms:
            hatoms.append(hetero)
            non_redundant_hs.append(h_index)

    for h_index in non_redundant_hs:
        new_mol = clone(mol)
        new_mol.DeleteAtom(new_mol.GetAtom(h_index))
        new_mol.SetTotalSpinMultiplicity(2)
        new_mol.SetTotalCharge(0)
        radicals.append(new_mol)
    return radicals

def neighbor(mol, index, hetero_indexs):
    q_atom = mol.GetAtom(index)
    for i in hetero_indexs:
        if mol.GetBond(index, i):
            #print "#%d is attached to atom #%d" % (index, i)
            return i

if __name__ == '__main__':
    input = sys.argv[1]

    obc = ob.OBConversion()
    obc.SetInAndOutFormats("smi", "mol")

    mol = ob.OBMol()
    next = obc.ReadString(mol, input)

    radicals = abs_h(mol)
    for r in radicals:
        print obc.WriteString(r)

GAMESSラッパー

Sat, 06 Nov 2010 18:57:39 +0900

量子化学計算できないケミストは、BLAST検索できない自称分子生物学者みたいなもんじゃねーの？

と常々思っている。

で、そういうヒトってなにに基づいてケミストリーしてんだろ？という素朴なギモンがあるのだが、ここではそういう部分は抜きにして、とりあえずBLASTのようにweb上から簡単に計算できるような量子化学計算サービスを用意すれば、リテラシーの問題で量子化学計算に手を出せなかった層に訴求できるんじゃないかな？といつも思っていたのだけど、所詮他人事なのでスルーしてきた。

でも、いろいろスルーできない状況になってきたので、余暇を利用して少しそこら辺も考えてみたというか、ハッカソンのテーマとして自分に課したら、結構はかどったので良かった、素晴らしい。

というわけで、ケモインフォクックブックを更新した。

次回はIzu.R #1として、発表付きのハッカソン形式としてやれればいいかなと思っている。

openbabel+GAMESSで計算できるようにする

Sat, 30 Oct 2010 21:35:05 +0900

Izu.R #0お疲れ様でした。いつもは22時前には寝るのだけど、2時半くらいまで起きてコード書いてた。たのしす。

今回が初めてで、Rらしく発表+LT的なほうがいいのか、ハッカソン的なもののほうがいいのか、わからなかったうえに初めての宿という不安はあったが、とりあえず集まってみて様子見しようという感じでしたが、蓋を開けてみれば、楽しい会になってくれて良かった。

次回もやる方向で考えているので、興味があれば。

今回のキーワード

Galaxy

pythonで書かれたワークフロー管理ツール。@bonohu さんに初期の導入のあたりを教えてもらいながら。みんな、ほとんどこれいじってたような。

Galaxy ツールを作る

ちなみにmacportsのpythonだと依存関係トラブリまくってうまく動かないので、潔く/usr/binのを使えという結論に落ち着いた。

これはもう少しいじる。次回にはなんか成果物をもっていきたいかな。

別トラックで流れていたセッションの情報交換

これはよい。あとで内容をきちんと把握する。あの、ユーザー会は一定の確率で面白い内容出てくるからよいです。面白さを面白いと感じるためには、面白さに対して感度をあげる努力を怠るとダメですな。

patent analysis(LDA) and spotfire

View more presentations from ybenjo.

ジャパネットムスカ

きっかけは覚えてないが、MADの方向に行ってきた。

ついでに、一人ふぁぼ界に引きずり込まれてた。

Python

Dive into Python

みんなのPythonは良書ですよね。

みんなのPython 改訂版
柴田淳
ソフトバンククリエイティブ / ￥ 2,940 ()
在庫あり。

そういえば、前日の二次会でもPython利用者率を上げていこうという結論で落ち着いた(Pythonistaの中で)

書いたもの

openbabelで荒く構造立ち上げて、その後量子化学計算で精密化して、Mulliken Chargeとか楽に求めたいこと多いので、サクサク出来るようにラッパー欲しかったので書いてた。ファイルの読み込みとかよくわからんバグが出ててはまったが、なんとか動くとこまでいった。

あとはもう少しちゃんと書いてクックブックにでもあげとく。

import openbabel as ob
from tempfile import mkstemp, mkdtemp
from os import removedirs, unlink, system, environ
import re
import os.path
import string
from random import choice

def randstr(n):
    "ランダムなファイル名を生成するため"
    return u''.join(choice('abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ') for i in xrange(n))

class Gamess(object):

    def __init__(self):
        self.tempdir = mkdtemp()
        self.gamess  = "/Users/kzfm/gamess/rungms"
        self.jobname = ''
        self.cwd    = os.getcwd()

    def calc(self, mol):
        self.jobname = randstr(6)
        obc = ob.OBConversion()
        obc.SetInFormat("gamout")

        gamin = self.write_file(mol)
        gamout = self.tempdir + "/" + self.jobname + ".out"
        os.chdir(self.tempdir)
        os.system("%s %s> %s  2> /dev/null" % (self.gamess, self.jobname, gamout))        

        # エラーが出たのでstringを渡したらなおった
        # TypeError: in method 'OBConversion_ReadFile', argument 3 of type 'std::string'
        new_mol = ob.OBMol()
        s = open(gamout).read()
        next = obc.ReadString(new_mol, s)

        os.chdir(self.cwd)
        unlink(gamin)
        unlink(gamout)
        return new_mol

    def header(self):
        h = """ $BASIS GBASIS=STO NGAUSS=3 $END
 $CONTRL SCFTYP=RHF RUNTYP=OPTIMIZE $END
 $STATPT OPTTOL=0.0001 NSTEP=20 $END"""

        return h

    def write_file(self,mol):
        obc = ob.OBConversion()
        obc.SetOutFormat("gamin")

        gamess_input_file = self.tempdir + "/" + self.jobname +".inp"
        gamess_input_str = obc.WriteString(mol)

        h = self.header()
        ng = gamess_input_str.replace(" $CONTRL COORD=CART UNITS=ANGS $END",h)

        with open(gamess_input_file, "w") as f:
            f.write(ng)
        return gamess_input_file

    def __del__(self):
        #print "remove " + self.tempdir
        removedirs(self.tempdir)

if __name__ == '__main__':

    g = Gamess()

    obc = ob.OBConversion()
    obc.SetInFormat("mol")

    mol = ob.OBMol()
    next = obc.ReadFile(mol,"test.mol")

    newmol = g.calc(mol)
    print [obatom.GetPartialCharge() for obatom in ob.OBMolAtomIter(newmol)]

GAMESSに関しては動く環境をもっているのが前提なので、そこをやりたい場合はこの本が良いです。

すぐできる量子化学計算ビギナーズマニュアル (KS化学専門書)

講談社 / ￥ 3,360 ()
在庫あり。

Pythonで演算子のオーバーライド

Thu, 05 Aug 2010 22:35:15 +0900

openbabelのPythonのとこ見てたらFingerprintsのとこで

>>> print fps[0] | fps[1] # Print the Tanimoto coefficient
0.3333

タニモト距離を求めるように、メソッドオーバーライドしてんのかと。コードを眺めた。

class Fingerprint(object):
    def __init__(self, fingerprint):
        self.fp = fingerprint
    def __or__(self, other):
        return ob.OBFingerprint.Tanimoto(self.fp, other.fp)
    @property 
    def bits(self):
            return _findbits(self.fp, ob.OBFingerprint.Getbitsperint())
    def __str__(self):
        return ", ".join([str(x) for x in self.fp])

参考

3.4.7 数値型をエミュレーションする

Python クックブック第2版
Alex Martelli,Anna Martelli Ravenscroft,David Ascher
オライリー・ジャパン / ￥ 4,410 ()
在庫あり。

openbabelのAddHydrogens

Tue, 12 Jan 2010 10:49:28 +0900

XYZでreadした分子(CC)もAddHydrogensできない

>>> from openbabel import *
>>> conv = OBConversion()
>>> conv.SetInFormat("xyz")
True
>>> mol = OBMol()
>>> conv.ReadFile(mol, "test.xyz")
True
>>> mol.NumAtoms()
2  
>>> mol.AddHydrogens()
True
>>> mol.NumAtoms()
2

でも分子中のそれぞれの原子を指定して

mol.AddHydrogens(atom)

としてやればきちんと水素が付加されることがわかったので、とりあえずスクリプト中でループをまわせばいいかな

openbabelで分子を構築するサンプル

Tue, 12 Jan 2010 10:49:15 +0900

pythonラッパーで。openbabelで気をつけないといけないところはAtomのインデックスは1から始まるのに対しBondのインデックスは0からはじまる。

from openbabel import *

mol = OBMol()

a = mol.NewAtom()
a.SetAtomicNum(6)
a.SetVector(0.0, 1.0, 2.0)

b = mol.NewAtom()
b.SetAtomicNum(6)
b.SetVector(1.0, 2.0, 2.0)

mol.AddBond(1, 2, 2) 
mol.SetTitle("Ethene")

obConversion = OBConversion()
obConversion.SetOutFormat("smi")

smiles_string = obConversion.WriteString(mol)

print smiles_string

あとAddHydrogensで水素を付与した場合の水素原子はSetVectorメソッドで座標が貼付けられないの。よくわからん。サンプルだときちんと水素がついてるので、なんかの属性をセットしとく必要があるのではなかろうか。

SDFからGASTON用のファイルを作る

Tue, 12 Jan 2010 10:43:01 +0900

gSpanやGASTONの入力はラベル付きのvertexとedgeなので、openbabelでsdfを読み込んで、GASTONのinputに変換するものを作ってみた。

import openbabel as ob

def convert(sdf):
    obc = ob.OBConversion()
    obc.SetInAndOutFormats('sdf','smi')

    mol = ob.OBMol()
    next = obc.ReadFile(mol,sdf)
    molnum = 0
    while next:
        # mol
        print "t # %d" % molnum
        # atom
        for i,atom in enumerate(ob.OBMolAtomIter(mol)): 
            print "v %d %d " % (i,atom.GetAtomicNum())

        # bond
        for i,bond in enumerate(ob.OBMolBondIter(mol)): 
            print "e %d %d %d" % (bond.GetBeginAtomIdx()-1,bond.GetEndAtomIdx()-1,bond.GetBondOrder())

        mol = ob.OBMol()
        next = obc.Read(mol)

        molnum += 1
    return True

if __name__ == "__main__":
    sdffile = 'pubchem_sample.sdf'

    convert(sdffile)

SDFはpubchemから適当に選んだがCID: 16757835は外しといた。

./gaston 11 pubchem_data pubchem_out

GASTONを実行した出力の一部

アロマティックなボンドは3みたいに別のラベルを与える必要がある気はするが、、、

後は、GASTONの出力ファイルから構造を構築するものを用意すれば、頻出する部分構造を扱えるようになる。

Molblasterを実装してみた

Tue, 12 Jan 2010 10:38:26 +0900

Molblasterっていうのは化学構造をランダムに切断していくっていう試行を何度か繰り返す方法で、その後フラグメントの頻度をシャノンエントロピーで評価したりするらしい。

シャノンエントロピーのほうに興味があったのとちょっと簡単に実装できそうだったのでpython+openbabelでやってみた。

構造をランダムに切断したい

まぁでも、構造活性相関にもっていくんだったら普通はベイズの方にいくよなぁ。多分著者の論文あると思うんだけど調べてない。