2012/01/28 09:10:47
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
去年書いた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を読んでます。二周目か三周目かわからんけど、何回読んでもこの本は楽しい。
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
量子化学計算できないケミストは、BLAST検索できない自称分子生物学者みたいなもんじゃねーの？
と常々思っている。
で、そういうヒトってなにに基づいてケミストリーしてんだろ？という素朴なギモンがあるのだが、ここではそういう部分は抜きにして、とりあえずBLASTのようにweb上から簡単に計算できるような量子化学計算サービスを用意すれば、リテラシーの問題で量子化学計算に手を出せなかった層に訴求できるんじゃないかな？といつも思っていたのだけど、所詮他人事なのでスルーしてきた。
でも、いろいろスルーできない状況になってきたので、余暇を利用して少しそこら辺も考えてみたというか、ハッカソン のテーマとして自分に課したら、結構はかどったので良かった、素晴らしい。
というわけで、ケモインフォクックブックを更新した 。
次回はIzu.R #1として、発表付きのハッカソン形式 としてやれればいいかなと思っている。
2010/10/30 21:35:05
Izu.R #0 お疲れ様でした。いつもは22時前には寝るのだけど、2時半くらいまで起きてコード書いてた。たのしす。
今回が初めてで、Rらしく発表+LT的なほうがいいのか、ハッカソン的なもののほうがいいのか、わからなかったうえに初めての宿という不安はあったが、とりあえず集まってみて様子見しようという感じでしたが、蓋を開けてみれば、楽しい会になってくれて良かった。
次回もやる方向で考えているので、興味があれば。
今回のキーワード
Galaxy
pythonで書かれたワークフロー管理ツール。@bonohu さん に初期の導入のあたりを教えてもらいながら。みんな、ほとんどこれいじってたような。
ちなみにmacportsのpythonだと依存関係トラブリまくってうまく動かないので、潔く/usr/binのを使えという結論に落ち着いた。
これはもう少しいじる。次回にはなんか成果物をもっていきたいかな。
別トラックで流れていたセッションの情報交換
これはよい。あとで内容をきちんと把握する。あの、ユーザー会は一定の確率で面白い内容出てくるからよいです。面白さを面白いと感じるためには、面白さに対して感度をあげる努力を怠るとダメですな。
ジャパネットムスカ
きっかけは覚えてないが、MADの方向に行ってきた。
ついでに、一人ふぁぼ界に引きずり込まれてた。
Python
みんなのPythonは良書ですよね。
そういえば、前日の二次会でも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
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])
参考
Python クックブック 第2版
Alex Martelli,Anna Martelli Ravenscroft,David Ascher
オライリー・ジャパン / ￥ 4,410 ()
在庫あり。
2010/01/12 10:49:28
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
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
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の出力ファイルから構造を構築するものを用意すれば、頻出する部分構造を扱えるようになる。
2010/01/12 10:38:26
Molblasterっていうのは化学構造をランダムに切断していくっていう試行を何度か繰り返す方法で、その後フラグメントの頻度をシャノンエントロピーで評価したりするらしい。
シャノンエントロピーのほうに興味があったのとちょっと簡単に実装できそうだったのでpython+openbabelでやってみた。
構造をランダムに切断したい
まぁでも、構造活性相関にもっていくんだったら普通はベイズの方にいくよなぁ。多分著者の論文あると思うんだけど調べてない。
エキスパートPythonプログラミング
みんなのPython 改訂版
すぐできる 量子化学計算ビギナーズマニュアル (KS化学専門書)
Python クックブック 第2版