24072015 三島 chemoinformatics q-chem
15092012 chemoinformatics q-chem
懇親会でpygamessじゃなくてもFacioで良かろうと言われたのだけど、対話環境よりも、スクリプトでやれたほうがいい場面が多いので書いておく。ちなみに僕はAvogadroユーザーです。
適当なsdfが手元になかったので
C CC CCC
というSMILESのリストを
babel -ismi carbons.smi -osdf carbons.sdf
とやってsdfにしたものを使いました。
sdfの分子群に対して計算したい場合、対話環境だとチマチマと作業しないといけないが、スクリプトだとforループを回せばよいし、結果がmoleculeオブジェクトで返ってくるので、そのあとのケモインフォマティクス的な取り回しが楽です。Gamessのアウトプット睨みながらコピペとか苦行だし、アウトプットの処理用のパーサー書くならinputから全てスクリプトで処理したほうがいいんじゃないかなと。
>>> import pybel >>> from pygamess import Gamess >>> g = Gamess() >>> mols = pybel.readfile('sdf', 'carbons.sdf') >>> for mol in mols: ... nmol = g.run(mol) ... print nmol, nmol.energy ... C -39.7265813363 CC -78.3052918313 CCC -116.885136991
ケミストのように化合物を個別のものとして一つ一つ丁寧に対応していくのか、情報論的にまとめて取り扱いたいのかとかそういうあたりで使い分ければいいのだと思う。オービタル見たい場合にはGUI必須だけど、エネルギーみたいな数字だけに関心があるのであればGUIは要らないかな。
バグフィックスしたり基底関数まわりを追加してバージョンをちょっとあげた。
一番大きな変更はデフォルトでrungmsを使わないようにしたことで 直接ddikickを叩くことにしたことか。それでgamessのパスをどうしようかなぁと。
とりあえずコンストラクタで直接指定してない場合は
それとは別に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みたいなノリで計算できるようにしたい。
pygamessがrungmsを設定してない場合でもマトモなエラーを吐かないとレスポンスを貰ったので、きちんと例外投げるように変更しておいた。ついでにドキュメントを追加してバーションを上げておいた。
pygamessはもともと、gamessのインプット作るのめんどくさいっていうのがモジュール作成のモチベーションだったので、gamess実行環境(rungms)が既にあるという前提だったが、最近pybelに対応したことにより、分子設計的な側面が強くなってしまったので、gamess実行環境もまとめて面倒みたほうがいいかなぁと。
rungmsは単なるシェルスクリプトなんだけど、ファイルを直接編集しないといけないので、動くようにするのがめんどくさい(色々なOSに対応するようにしているのでごちゃごちゃしている)。普通に動かせる最小のスクリプトはどんな感じかなぁと1000行超えのスクリプトを削って行ったら200行くらいになったが、残ったコードのほとんどがsetenvだった。どんだけ環境変数好きやねん?と。
これをpythonで書きなおしたら60行くらい。
import os import sys import socket from shutil import copyfile, rmtree from tempfile import mkstemp, mkdtemp scr = mkdtemp() job = sys.argv[1] gamess_path = "/usr/local/gamess" ddikick = os.path.join(gamess_path, "ddikick.x") gamess = os.path.join(gamess_path, "gamess.Jan122009R1.x") hostname = socket.gethostname() setenv_data = [ (" MAKEFP", "efp"), ("GAMMA", "gamma"), ("TRAJECT", "trj"), ("RESTART", "rst"), (" PUNCH", "dat"), (" INPUT", "F05"), (" AOINTS", "F08"), (" MOINTS", "F09"), ("DICTNRY", "F10"), ("DRTFILE", "F11"), ("CIVECTR", "F12"), ("CASINTS", "F13"), (" CIINTS", "F14"), (" WORK15", "F15"), (" WORK16", "F16"), ("CSFSAVE", "F17"), ("FOCKDER", "F18"), (" WORK19", "F19"), (" DASORT", "F20"), ("DFTINTS", "F21"), ("DFTGRID", "F22"), (" JKFILE", "F23"), (" ORDINT", "F24"), (" EFPIND", "F25"), ("SVPWRK1", "F26"), ("SVPWRK2", "F27"), (" MLTPL", "F28"), (" MLTPLT", "F29"), (" DAFL30", "F30"), (" SOINTX", "F31"), (" SOINTY", "F32"), (" SOINTZ", "F33"), (" SORESC", "F34"), ("GCILIST", "F37"), ("HESSIAN", "F38"), ("QMMMTEI", "F39"), ("SOCCDAT", "F40"), (" AABB41", "F41"), (" BBAA42", "F42"), (" BBBB43", "F43"), (" MCQD50", "F50"), (" MCQD51", "F51"), (" MCQD52", "F52"), (" MCQD53", "F53"), (" MCQD54", "F54"), (" MCQD55", "F55"), (" MCQD56", "F56"), (" MCQD57", "F57"), (" MCQD58", "F58"), (" MCQD59", "F59"), (" MCQD60", "F60"), ("NMRINT1", "F61"), ("NMRINT2", "F62"), ("NMRINT3", "F63"), ("NMRINT4", "F64"), ("NMRINT5", "F65"), ("NMRINT6", "F66"), ("ELNUINT", "F67"), ("NUNUINT", "F68"), (" NUMOIN", "F69"), (" GMCREF", "F70"), (" GMCO2R", "F71"), (" GMCROC", "F72"), (" GMCOOC", "F73"), (" GMCCC0", "F74"), (" GMCHMA", "F75"), (" GMCEI1", "F76"), (" GMCEI2", "F77"), (" GMCEOB", "F78"), (" GMCEDT", "F79"), (" GMCERF", "F80"), (" GMCHCR", "F81"), (" GMCGJK", "F82"), (" GMCGAI", "F83"), (" GMCGEO", "F84"), (" GMCTE1", "F85"), (" GMCTE2", "F86"), (" GMCHEF", "F87"), (" GMCMOL", "F88"), (" GMCMOS", "F89"), (" GMCWGT", "F90"), (" GMCRM2", "F91"), (" GMCRM1", "F92"), (" GMCR00", "F93"), (" GMCRP1", "F94"), (" GMCRP2", "F95"), (" GMCVEF", "F96"), (" GMCDIN", "F97"), (" GMC2SZ", "F98"), (" GMCCCS", "F99") ] for e in setenv_data: os.environ[e[0].strip()] = "%s/%s.%s" %(scr, job, e[1]) os.environ["ERICFMT"] = os.path.join(gamess_path, "ericfmt.dat") os.environ["MCPPATH"] = os.path.join(gamess_path, "mcpdata") os.environ["EXTBAS"] = "/dev/null" os.environ["NUCBAS"] = "/dev/null" src = job + ".inp" dest = os.path.join(scr,job) + ".F05" copyfile(src,dest) exec_string = "%s %s %s -ddi 1 1 %s -scr %s > t.out" % (ddikick, gamess, job, hostname, scr) os.system(exec_string) rmtree(scr)
結局ユーザーが指定しないといけない変数ってGamessのpathくらいだった。
17032012 chemoinformatics q-chem Python
昨日は製薬業界の集まりがあって、他社のヒトと少し話す機会があって、LBDDどういう感じですかねと言われて、量子化学計算に真面目に取り組んでますよーっていう話から、ケミストは電子吸引基とか供与基とか言う割に計算して確かめようとしないんですよねー(そうですよねー)っていう流れになったので、そちらのケミストは計算するんですかねー?って聞いたら、するヒトはするし、しないヒトはしないっていう答えが返ってきた。想定通り。
もう少し知りたいのは、ちゃんと計算するケミストは、結局ただの計算マニアで結局普通のヒトなのか、それとも論理的で優秀な傾向が強いのか?ということかな。
さて、pygamessをpybelに対応させたので、より簡潔に書けるようになった。
>>> from pygamess import Gamess >>> import pybel >>> g = Gamess() >>> g.run_type('optimize') >>> mol = pybel.readstring('smi','C') >>> mol.make3D() >>> optimized_mol = g.run(mol) >>> optimized_mol.energy -37.0895866208
やっぱコンストラクタに引数をわたせるようにするべきだよなぁ。
量子化学計算したいならこれかな。
07032012 chemoinformatics q-chem Python
Gamessのpythonラッパーであるpygamessを0.2.0にアップデートした。
pybelも使えるようになったので、よりpythonicにコードを書けるようになっていい感じ。 構造さえ用意すれば、こんなに簡単に量子化学計算ができて、量子化学計算にありがちなごちゃごちゃしたインプットファイルの作成から解放される。
>>> import pybel, pygamess >>> g = pygamess.Gamess() >>> mol = pybel.readfile("mol", "examples/ethane.mol").next() >>> nmol = g.run(mol) >>> nmol.energy -78.30530748
それから、テストをpyvowsに変えてから快適です。自分にはBDDはあっているなぁ。
ドキュメントをきちんと書く
30112011 q-chem
EC2でFMO計算できたら、超速で楽しいんだろうなと、ずっと思っていたが、この前の勉強会でイメージが湧いたのでAmazon EC2の価格表見ていたら超やばい。
物理的なハードウェアがなくても、安く早く計算できそうですよ。ハイメモリ クアドラプル エクストララージ インスタンスを4つくらい使えば大抵のタンパク質は数時間で計算終わりそうな気がするなぁ。
というわけでFMO計算とSBDDとEC2というノウハウを組み合わせればニッチな受託モデルができそうな気が。
でも、量子化学計算を活用したStructure Based Drug Designってニーズあんのかなぁ?ウェットのケミストの理解を得られない気もするが。
13072011 chemoinformatics q-chem Python
AM1とかちょっと使いたかったので使えるようにしたけど、そうするとエネルギーが取り込めなくて困る。
あとイオン化ポテンシャルも取り込みたいのでそこら辺も対応する予定。
論文読んでたらハモンド仮説ってのがあって、あれってなんだったかなとググッたら自分の持ってる本が出てきたので、家に帰ってから読んでみたら、かなり面白かったのでこの本はオススメ。
メカニズムベースで創薬研究(特に合成)やろうと思ったら量子化学の基本的なところは外せないと思うんだよね。KKDってのはそれを超えたところにあるわけであって、宝くじ買って当てようとすることをセレンディピティとは言わないでしょう。
去年からダラダラと作っていたPythonのGAMESSラッパーをPyPIにあげてみたけど、初めてなので、なんかおかしい部分があったら指摘してもらえるとありがたいです。
登録は以下のサイトを参考にした
実際にやってみて分かりづらかったところはやっぱライセンスかなぁ。逆にclassifierのところはlist見ながら該当する部分をコピペしてけばいいだけだったので思っていたよりは楽だったな。
ROHFで。
! File created by the GAMESS Input Deck Generator Plugin for Avogadro $BASIS GBASIS=STO NGAUSS=3 $END $CONTRL SCFTYP=ROHF RUNTYP=OPTIMIZE MULT=2 $END $STATPT OPTTOL=0.0001 NSTEP=20 $END $DATA Title C1 C 6.0 -1.60293 -0.39674 -0.00000 H 1.0 -0.53293 -0.39674 -0.00000 H 1.0 -1.97337 -0.92091 0.85611 H 1.0 -1.95466 -0.89443 -0.87948 $END
平面になった。
TODO: pythonのGAMESSラッパーをPyPIに登録する