Drkcore

高速文字列解析の流れで、タイミングよく統数研チャンネルの資料が公開されていて、最後のほうのスライド(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 )

10012013 chemoinformatics Python openbabel

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

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

31032012 q-chem Python openbabel

pygamess 0.3.0

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

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みたいなノリで計算できるようにしたい。

30032012 Python openbabel

素晴らしいスピン多重度

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

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

29032012 Python openbabel

pybelでanilineをnitrenium ionに変換する

なぜ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のオブジェクトを直接触る必要があるのであまりスマートではない気がする。

23022012 chemoinformatics Python openbabel

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

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パラメータを使って最急降下法で構造最適化をしている。

28012012 chemoinformatics Python openbabel

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

22062011 chemoinformatics Python openbabel GAMESS

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

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

18062011 chemoinformatics openbabel

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

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

06112010 chemoinformatics openbabel GAMESS

GAMESSラッパー

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

と常々思っている。

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

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

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

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

もう5年目(wishlistありマス♡)
最近はPythonとDeepLearning
日本酒自粛中
ドラムンベースからミニマルまで
ポケモンGOゆるめ