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

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

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

今回のキーワード

Galaxy

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

• Galaxy ツールを作る

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

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

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

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

ジャパネットムスカ

きっかけは覚えてないが、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 ()
在庫あり。

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 ()
在庫あり。

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)

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

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

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を実行した出力の一部

# 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の出力ファイルから構造を構築するものを用意すれば、頻出する部分構造を扱えるようになる。

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

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

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

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

openbabelのringsearchがちょっと使いにくい気がする。部分構造をmolオブジェクトで返すメソッドがあっても良いんではないか？

from openbabel import *

obconversion = OBConversion()
obconversion.SetInFormat("smi")
obconversion.SetOutFormat("smi")
obmol = OBMol()

notatend = obconversion.ReadString(obmol,"c1ccccc1CC(=O)C2CCC2C(=O)O")

for r in obmol.GetSSSR():
    for a in  r._path:
        print a

これだと原子のIDしか返さない。

でも縮環チェックしたいんだったら、原子のリストがかぶるかどうかをリングの構成原子のIDをつかって調べればいいのでminimumの環からmaximumの環にするのは難しくはないか。

アトムのリストを渡したら部分構造検索を実行して、結果としてマッチした構造をmolとして返すメソッドがあればいいのか。

Open Babel 2.2.0にバージョンがあがったのでインストールをした。ソースから。

同時に、perl,python,rubyのバインディングもコンパイルしてインストールしておく。

WindowsにおいてあるPlutoのファイル群を持ってくる

Windowsで開発しているPlutoをmacでもいじれるようにする。windowsのほうはhg serveと打てばhttpサーバーが立ち上がり、port8000番でアクセスできるようになる。

macのほうはディレクトリを作って初期化してpullしてupdate

mkdir Pluto
cd Pluto
hg init
hg pull http://192.168.XXX.XXX:8000/
hg update

これでOK

openbabelの分子重ね合わせツールのobfitをperlのモジュールにできんもんかな？とソース(C++)を眺めていたら、QTRFIT というアルゴリズム（実装はC++)を使っているようだ。

結局最大の重ねあわせを求めるのを固有値問題に落とし込んでるようだ。

QTRFIT

The eigenvector associated with the largest eigenvalue gives the desired quaternion. To generate a rotation matrix we consider the effects of the quaternion on the cartesian unit vectors

で、それをJacobi法で解いていると。感覚的にはそうなんだろなぁと思うが、ここら辺の理解度はちょっとあやしい。