openbabel / Drkcore

Graph-indexing wavelet tree (gWT)のグラフ検索が超速い

Fri, 11 Jan 2013 05:54:20 +0919

高速文字列解析の流れで、タイミングよく統数研チャンネルの資料が公開されていて、最後のほうのスライド(65枚目)にグラフの類似度検索の実装としてgwtが紹介されていたので遊んでみた。

データベースは、ChEMBLから取ってきた。なぜpubchemじゃないかというと面倒なので略。

入力のフォーマットがgSpanとかGASTONと一緒だったので、昔書いたのを流用しようとしたんだが、昔書いたpythonスクリプトはどうもOBMolAtomIterあたりがバグっていて動かないようだ(OBMolBondIterも動かなかった。)

File "sdf2gsp", line 17, in convert
  for i,atom in enumerate(ob.OBMolAtomIter(mol)): 
File "/usr/local/Cellar/python/2.7.3/lib/python2.7/site-packages/openbabel.py", line 4499, in __init__
  if not self.iter.__bool__():
AttributeError: '_OBMolAtomIter' object has no attribute '__bool__'

しょうがないのでiterateしない版で書き換えた。

#!/usr/bin/env python
import sys
import pybel
import openbabel as ob

def convert(sdf):
    mols = pybel.readfile('sdf',sdf)
    for n, mol in enumerate(mols, 1):
        print "t # %d" % n
        for i in range(mol.OBMol.NumAtoms()):
            print "v %d %d " % (i+1, mol.OBMol.GetAtom(i+1).GetAtomicNum())
        for j in range(mol.OBMol.NumBonds()):
            b = mol.OBMol.GetBond(j)
            print "e %d %d %d" % (b.GetBeginAtomIdx(),
                                  b.GetEndAtomIdx(),
                                  b.GetBondOrder())

if __name__ == "__main__":
    sdffile = sys.argv[1]
    convert(sdffile)

ちなみにこの形式に名前ついているんだろうか？babelとかでサポートされないと使い勝手が悪いんだけどなぁ。

これで、ChEMBLの137万件をgspにコンバートしてみたけど遅すぎるので後でHaskellで書きなおす。sdf単に舐めるだけだからOBMolに変換する必要もないし。

できたchembl_14.gspからindexファイルをつくる。600M(137万化合物)で10分ぐらいかかった@MBA。

./gwt-build -iteration 2 chembl_14.gsp chembl14

検索はみんな大好きImatinibをgspにしたものを使った。

$ ./gwt-search -kthreshold 0.8 chembl14 query.gsp
reading index file:chembl14
end reading
start building rank dictionary
end building rank dictionary
start transferring labels
end transferring labels
id:0 1080597:0.809449 712949:0.838758 413038:0.847826 \
394513:0.826087 391369:0.819028 179393:0.815124 \
175534:0.899647 48523:0.804348 47164:1.000000 16405:0.802727 
cpu time (sec):0.712119
mean cpu time (sec):0.712119

と5,6秒程度で計算終了。ウリィィィィって言いたくなるほど圧倒的に速い。

bayonも感激したけど、それ以上にgwtはすごいなぁと思った。

高速文字列解析の世界――データ圧縮・全文検索・テキストマイニング (確率と情報の科学)
岡野原大輔
岩波書店 / 3150円 ( 2012-12-27 )

OS10.8.2にopenbabel2.3.2をインストール

Thu, 10 Jan 2013 20:06:02 +0919

2.3.1のままだったので、昔のログを見ながらインストール。

が、古いswig(/usr/binにあるやつ)を見ていてエラーになるので、SWIGの場所をcmakeに知らせたいのだがオプションがわからない。

調べたトコロ

cmake -L

で一覧が出るらしい。必要なのはSWIG_EXECUTABLEなので

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 \
-DSWIG_EXECUTABLE=/usr/loca/bin/swig

でOK

pygamess 0.3.0

Sat, 31 Mar 2012 13:26:32 +0919

バグフィックスしたり基底関数まわりを追加してバージョンをちょっとあげた。

pygamess

一番大きな変更はデフォルトでrungmsを使わないようにしたことで直接ddikickを叩くことにしたことか。それでgamessのパスをどうしようかなぁと。

とりあえずコンストラクタで直接指定してない場合は

環境変数GAMESS_PATHを調べる
PATHの中からddikick.xのあるディレクトリを探して、それを選択する

それとは別にrungmsのpathを探したりもしている(僕の場合はrungmsを/usr/local/binに置いている)のでちょっとくどいやり方になってしまった。

それから、コンストラクタにオプションを渡せるようにしたんだけど、これだと、gamessのインプットいじるのと変わらないので、もう少しスマートな指定方法も用意したほうがいいかなぁと思った。

水分子を基底関数を変えながらエネルギーを計算する例。このくらいだったら対話環境でサクッとうごく。

>>> import pybel
>>> from pygamess import Gamess
>>> g = Gamess()
>>> mol = pybel.readstring('smi', 'O')
>>> mol.make3D()
>>> g.run(mol).energy
-74.96450135
>>> g.run_type('optimize')
>>> g.run(mol).energy
-74.9659012146
>>> g.basis_set('3-21G')
>>> g.run(mol).energy
-75.585959758
>>> g.basis_set('6-31G')
>>> g.run(mol).energy
-75.9853591564
>>> g.basis_set('6-311G')
>>> g.run(mol).energy
-76.0109546389
>>> g.basis_set('6-31G*')
>>> g.run(mol).energy
-76.0107465155
>>> g.basis_set('6-31G**')
>>> g.run(mol).energy
-76.0236150193

最近ちょっと困っているのが、構造最適化で収束しないのをどうやって収束させるかという。6-31G*くらいでmopacみたいなノリで計算できるようにしたい。

素晴らしいスピン多重度

Fri, 30 Mar 2012 19:25:04 +0919

openbabelのコンバーターがバグってる。

>>> import pybel
>>> from pygamess import Gamess
>>> g = Gamess()
>>> mol = pybel.readstring('smi','C')
>>> mol.make3D()
>>> nmol = g.run(mol)
>>> nmol.OBMol.GetTotalSpinMultiplicity()
1882919584

pybelでanilineをnitrenium ionに変換する

Thu, 29 Mar 2012 19:08:54 +0919

なぜnitrenium ionにする必要があるのかはここではふれないので、ココらへんを見てください。

import pybel
smarts = pybel.Smarts("c[#7]([#1])[#1]")

mol = pybel.readstring('smi', 'c1ccccc1N')
mol.make3D()
matches = smarts.findall(mol)
mol.OBMol.DeleteAtom(mol.atoms[matches[0][2]-1].OBAtom)
mol.OBMol.SetTotalCharge(1)
print mol.write('mol')

分子に修飾をしたい時はpybelにメソッドがないので、OpenBabelのオブジェクトを直接触る必要があるのであまりスマートではない気がする。

pybelで立体構造を立ち上げる

Thu, 23 Feb 2012 20:11:11 +0919

pybel(openbabel)で立体構造の立ち上げを行うのは簡単だ。

import pybel
mol = pybel.readstring("smi", "c1ccccc1C")
mol.make3D()
mol.write("mol")

コードを眺めていたら構造をたちあげるメソッドは

make3d
localopt
globalopt # ただし2.3.1ではコメントアウトされてた

make3dもglobaloptもlocaloptのstep数を変えて呼び出している。

localoptメソッドはデフォルトではmmff94パラメータを使って最急降下法で構造最適化をしている。

openbabel-2.3.1が出てますね。

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

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 +0919

去年書いた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 +0919

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

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

と常々思っている。

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

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

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

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