2011/09/23 18:24:28
光異性化とかAMES予測とか代謝予測とか、量子化学計算が創薬シーンで果たす役割が大きくなっているのは、探索創薬自体が発見学というよりはメカニズムベースでモノを考えるようになってきているからかなぁと。それからFMOなんかも重要な技術ですね。
さて、ちょっと励起状態での構造最適化計算が必要になったので、pygamess でCIS計算できるようにしておきました。
test用にGAMESSのEXAM34のホルムアルデヒドの励起状態計算のサンプルを使っています。GAMESSで計算する場合にはpc-chem.infoのホルムアルデヒドの励起状態計算 が参考になります。こういう良質なコンテンツもっと増えてくんないかな。
input用のmol
exam34_energy . log
OpenBabel09231115413D
4 3 0 0 0 0 0 0 0 0999 V2000
0.0100 - 0.8670 0.0000 O 0 0 0 0 0 0 0 0 0 0 0 0
0.0000 0.3455 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0
- 0.0100 0.9296 - 0.9377 H 0 0 0 0 0 0 0 0 0 0 0 0
- 0.0100 0.9296 0.9377 H 0 0 0 0 0 0 0 0 0 0 0 0
1 2 2 0 0 0 0
2 4 1 0 0 0 0
3 2 1 0 0 0 0
M END
計算用スクリプト
import pygamess
import openbabel as ob
g = pygamess . Gamess ()
obc = ob . OBConversion ()
obc . SetInAndOutFormats ( "mol" , "mol" )
mol = ob . OBMol ()
obc . ReadFile ( mol , "h2co.mol" )
g . contrl [ 'cityp' ] = 'cis'
g . run_type ( 'optimize' )
g . basis_type ( '631+gdp' )
g . gamess_input ( mol )
newmol = g . run ( mol )
print newmol . GetEnergy ()
obc . WriteFile ( newmol , 'h2co_singlet.mol' )
励起状態のホルムアルデヒドの安定構造は若干ピラミッド型の構造を取るって知ってた？
ところで、光異性化のサンプルとして面白いのはやっぱスチルベンかなぁと思ったんだけど、HOMO-LUMOの2電子励起を考慮しないといけないらしい ので、CISじゃ計算できないじゃんと。
もう少し面白いサンプルないかなぁ。
2011/08/06 16:37:39
生物学的等価体として知られているものにカルボン酸とテトラゾールがありますね。
医薬化学において、テトラゾール環はカルボン酸の等価体と見なされ、医薬品の部分構造に汎用されている。これは前述のようにpKaがほぼ等しいため、
Wikipedia
テトラゾールってこんな構造だからpKaは同じだとしてもアニオン化した状態で負電荷どこにたまるのかなという素朴な疑問を持っていた。
ま、計算すればいいだけなので、pygamess を使って計算してみた。以下コード。チャージがcontrlセクションに入力するっていうのが気持ち悪い。
あと、最適化がなかなか終わらないのでmaxitとnstepは増やした。デフォルトをこの値にしてもいいかもしれない。
from pygamess import Gamess , GamessError
import openbabel as ob
g = Gamess ()
g . contrl [ 'maxit' ] = 200
g . contrl [ 'icharg' ] = - 1
g . statpt [ 'nstep' ] = 300
g . system [ 'mwords' ] = 80
g . basis_type ( '631g' )
g . run_type ( 'optimize' )
obc = ob . OBConversion ()
obc . SetInAndOutFormats ( "mol" , "mol" )
mol = ob . OBMol ()
obc . ReadFile ( mol , "tetrazole.mol" )
print g . gamess_input ( mol )
try :
newmol = g . run ( mol )
except GamessError , gerr :
print gerr . value
else :
for obatom in ob . OBMolAtomIter ( newmol ):
print obatom . GetIdx (), obatom . GetType (), obatom . GetPartialCharge ()
いろいろ試していたので、GAMESSインプットを出力するようにしてある。
結果を見ると、確かに6,9番目の窒素のマイナスチャージが大きくなっているので、カルボン酸等価なんだろうなと納得したのであった。
$ /opt/local/bin/python optimize.py
$contrl runtyp = optimize scftyp = rhf icharg = -1 maxit = 200 mult = 1 $end
$basis gbasis = N31 ndfunc = 1 ngauss = 6 $end
$system mwords = 80 $end
$statpt opttol = 0.0001 nstep = 300 $end
$DATA
C1
C 6.0 -5.1500000000 0.5694000000 0.3351000000
C 6.0 -3.6632000000 0.3894000000 0.3044000000
H 1.0 -5.3940000000 1.6144000000 0.4937000000
H 1.0 -5.6014000000 -0.0166000000 1.1314000000
H 1.0 -5.6123000000 0.2533000000 -0.5982000000
N 7.0 -2.8049000000 1.3874000000 0.2522000000
N 7.0 -1.6217000000 0.7950000000 0.2303000000
N 7.0 -1.7700000000 -0.4762000000 0.2679000000
N 7.0 -3.0563000000 -0.7807000000 0.3162000000
$END
1 C3 -0.5149
2 Car 0.437823
3 HC 0.153658
4 HC 0.142616
5 HC 0.138266
6 Nar -0.487902
7 Nar -0.192393
8 Nar -0.188063
9 Nar -0.489106
量子化学計算っていうのは意外にトライアンドエラーが多い。だからといってGUIのコンソールからちまちまと作業すると時間が取られて困るし、作業ログが残りにくいので、こういったスクリプトで幾つかの条件を用意してforループでグルグルまわして計算させれば手順がきちんとドキュメントとして残るのでよい。
2011/06/25 10:11:48
去年からダラダラと作っていたPythonのGAMESSラッパー をPyPIにあげてみたけど、初めてなので、なんかおかしい部分があったら指摘してもらえるとありがたいです。
登録は以下のサイトを参考にした
実際にやってみて分かりづらかったところはやっぱライセンスかなぁ。逆にclassifierのところはlist 見ながら該当する部分をコピペしてけばいいだけだったので思っていたよりは楽だったな。
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 10:34:20
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に登録する
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に関しては動く環境をもっているのが前提なので、そこをやりたい場合はこの本が良いです。
2009/07/04 20:21:26
Cation-πってどのくらい安定化すんだろうか?と思っていたので計算してみる。
RCSBを適当にあさって見つけた、1R9L のBETとTRP188を切り出してモデル化する。
で、GAMESSで計算。
最初は6-31G+(d,p)で最適化しようと思ったのだけど、macbookだと全然終わらん。なのでSTO3Gでいいやと。TRPは無駄に計算量が増えるのだけどPHE,TYRあたりでモデル化できるのであればmacbookでも6-31G+(d,p)で計算できるような気がする。
------------------------------------
RESULTS OF KITAURA-MOROKUMA ANALYSIS
------------------------------------
HARTREE KCAL/MOLE
ELECTROSTATIC ENERGY ES= -0.011414 -7.16
EXCHANGE REPULSION ENERGY EX= 0.008788 5.51
POLARIZATION ENERGY PL= -0.001679 -1.05
CHARGE TRANSFER ENERGY CT= -0.007053 -4.43
HIGH ORDER COUPLING ENERGY MIX= 0.000166 0.10
TOTAL INTERACTION ENERGY, DELTA-E= -0.011193 -7.02
水素結合と同じくらいの強さかな。
以下、STO-3Gで構造最適化した座標を使ってつくったインプット
$BASIS GBASIS=STO NGAUSS=3 NDFUNC=1 $END
$CONTRL SCFTYP=RHF RUNTYP=MOROKUMA MAXIT=200 ICHARG=1 MULT=1 $END
$MOROKM IATM(1)=17 ICHM(1)=1 $END
$DATA
W188_BET.out
C1
N 7.0 26.56861 22.88472 72.47206
C 6.0 25.05591 23.05763 72.38527
C 6.0 27.18276 23.04159 71.08205
C 6.0 27.14794 23.94295 73.40648
C 6.0 26.89124 21.49601 73.01604
H 1.0 24.65761 22.29324 71.71973
H 1.0 24.84300 24.05067 71.99229
H 1.0 26.95089 24.04166 70.71972
H 1.0 26.73892 22.29197 70.42917
H 1.0 28.26387 22.89551 71.15590
H 1.0 26.89847 24.92550 73.00913
H 1.0 28.22884 23.81270 73.44995
H 1.0 26.70726 23.81001 74.39332
H 1.0 26.46855 20.75492 72.33939
H 1.0 26.44825 21.40230 74.00640
H 1.0 27.97398 21.38820 73.06968
H 1.0 24.63622 22.94695 73.38408
C 6.0 30.61324 20.49644 72.40614
C 6.0 30.47183 21.03406 71.19420
C 6.0 30.79531 21.59051 73.36029
C 6.0 30.75216 22.77934 72.62888
C 6.0 31.00908 21.63948 74.74382
N 7.0 30.43745 22.45766 71.27884
C 6.0 30.92671 24.02761 73.23103
C 6.0 31.18435 22.86580 75.34238
C 6.0 31.14962 24.05129 74.58925
H 1.0 30.33993 20.55500 70.23268
H 1.0 31.04481 20.72825 75.32842
H 1.0 30.86231 23.00634 70.52231
H 1.0 30.90815 24.93911 72.64846
H 1.0 31.36065 22.92555 76.40870
H 1.0 31.30378 24.99933 75.08996
H 1.0 30.61446 19.44515 72.64676
$END
2009/06/17 21:52:08
ありがちなエラー。
*** ERROR! THERE ARE NOT 5 OR 6 TRANS/ROT MODES NUM
そういう場合には
$statpt projct=.f. $end
と入れとく
2009/06/11 23:26:29
こんな感じ？
$BASIS GBASIS=N31 NGAUSS=6 NDFUNC=1 DIFFSP=.TRUE. $END
$CONTRL SCFTYP=RHF RUNTYP=OPTIMIZE DFTTYP=B3LYP NUMGRD=.T. $END
$SVP $END
$STATPT OPTTOL=0.0001 NSTEP=200 $END
$DATA
Title
C1
O 8.0 -4.21789 1.92792 0.00000
H 1.0 -3.25037 1.97435 0.00000
H 1.0 -4.49708 2.85544 0.00000
$END
多分だめだな。あとでやり直す
エキスパートPythonプログラミング
みんなのPython 改訂版
すぐできる 量子化学計算ビギナーズマニュアル (KS化学専門書)