drkcore

2012/01/28 09:10:47

openbabel-2.3.1が出てますね。

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

2011/06/22 19:07:02

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

去年書いた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を読んでます。二周目か三周目かわからんけど、何回読んでもこの本は楽しい。

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


2011/06/18 15:08:48

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

共有結合先の原子をチェックして冗長なものが出ないようにした。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)

2010/11/06 18:57:39

GAMESSラッパー

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

と常々思っている。

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

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

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

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

2010/10/30 21:35:05

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

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

1288429089 1288429091

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

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

今回のキーワード

Galaxy

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

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

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

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

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

ジャパネットムスカ

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

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

Python

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

ProductName みんなの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に関しては動く環境をもっているのが前提なので、そこをやりたい場合はこの本が良いです。

2010/08/05 22:35:15

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

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])

参考

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

2010/01/12 10:49:28

openbabelのAddHydrogens

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)

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

2010/01/12 10:49:15

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

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メソッドで座標が貼付けられないの。よくわからん。サンプルだときちんと水素がついてるので、なんかの属性をセットしとく必要があるのではなかろうか。

2010/01/12 10:43:01

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

gSpanGASTONの入力はラベル付きの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を実行した出力の一部

# 14
t 163
v 0 6
v 1 6
v 2 1
e 0 1 2
e 1 2 1
# 12
t 164
v 0 6
v 1 6
v 2 1
v 3 1
e 0 1 2
e 0 3 1
e 1 2 1

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

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

2010/01/12 10:38:26

Molblasterを実装してみた

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

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

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

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