光異性化とか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じゃ計算できないじゃんと。

もう少し面白いサンプルないかなぁ。

生物学的等価体として知られているものにカルボン酸とテトラゾールがありますね。

医薬化学において、テトラゾール環はカルボン酸の等価体と見なされ、医薬品の部分構造に汎用されている。これは前述のように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ループでグルグルまわして計算させれば手順がきちんとドキュメントとして残るのでよい。

去年からダラダラと作っていたPythonのGAMESSラッパーをPyPIにあげてみたけど、初めてなので、なんかおかしい部分があったら指摘してもらえるとありがたいです。

• pygamess 0.1

登録は以下のサイトを参考にした

• PyPIデビュー
• 簡単なPythonプロジェクトをPyPiに登録する

実際にやってみて分かりづらかったところはやっぱライセンスかなぁ。逆にclassifierのところはlist見ながら該当する部分をコピペしてけばいいだけだったので思っていたよりは楽だったな。

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

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

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に登録する

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

と常々思っている。

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

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

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

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

Izu.R #0お疲れ様でした。いつもは22時前には寝るのだけど、2時半くらいまで起きてコード書いてた。たのしす。

今回が初めてで、Rらしく発表+LT的なほうがいいのか、ハッカソン的なもののほうがいいのか、わからなかったうえに初めての宿という不安はあったが、とりあえず集まってみて様子見しようという感じでしたが、蓋を開けてみれば、楽しい会になってくれて良かった。

次回もやる方向で考えているので、興味があれば。

今回のキーワード

Galaxy

pythonで書かれたワークフロー管理ツール。@bonohu さんに初期の導入のあたりを教えてもらいながら。みんな、ほとんどこれいじってたような。

• Galaxy ツールを作る

ちなみにmacportsのpythonだと依存関係トラブリまくってうまく動かないので、潔く/usr/binのを使えという結論に落ち着いた。

これはもう少しいじる。次回にはなんか成果物をもっていきたいかな。

別トラックで流れていたセッションの情報交換

これはよい。あとで内容をきちんと把握する。あの、ユーザー会は一定の確率で面白い内容出てくるからよいです。面白さを面白いと感じるためには、面白さに対して感度をあげる努力を怠るとダメですな。

ジャパネットムスカ

きっかけは覚えてないが、MADの方向に行ってきた。

ついでに、一人ふぁぼ界に引きずり込まれてた。

Python

• Dive into Python

みんなのPythonは良書ですよね。

みんなのPython 改訂版
柴田 淳
ソフトバンククリエイティブ / ￥ 2,940 ()
在庫あり。

そういえば、前日の二次会でも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に関しては動く環境をもっているのが前提なので、そこをやりたい場合はこの本が良いです。

すぐできる 量子化学計算ビギナーズマニュアル (KS化学専門書)

講談社 / ￥ 3,360 ()
在庫あり。

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

すぐできる 量子化学計算ビギナーズマニュアル (KS化学専門書)

講談社 / ￥ 3,360 ()
在庫あり。

ありがちなエラー。

 *** ERROR! THERE ARE NOT 5 OR 6 TRANS/ROT MODES NUM 

そういう場合には

$statpt projct=.f. $end

と入れとく

こんな感じ？

 $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

多分だめだな。あとでやり直す

すぐできる 量子化学計算ビギナーズマニュアル (KS化学専門書)

講談社 / ￥ 3,360 ()
在庫あり。