Cation-π相互作用をSAPT0してみた。

別件で調べ物をしていたら10年くらい前に検討したエントリを見つけたので、これをpsi4のSAPTでやったらどういう結果になるかな?と。

今回psi4コマンドで流すのでinputはこんな感じです(いずれはpsikitでサクッと計算できるようにしたい)。EmacsやめちゃったからVimでgamessのインプットいじってつくった。繰り返しの技を覚えたけど、こういうの手になじまないとめんどくさいですね。

molecule cation_pi {
     1 1
     N 26.56861  22.88472  72.47206
     C 25.05591  23.05763  72.38527
     C 27.18276  23.04159  71.08205
     C 27.14794  23.94295  73.40648
     C 26.89124  21.49601  73.01604
     H 24.65761  22.29324  71.71973
     H 24.84300  24.05067  71.99229
     H 26.95089  24.04166  70.71972
     H 26.73892  22.29197  70.42917
     H 28.26387  22.89551  71.15590
     H 26.89847  24.92550  73.00913
     H 28.22884  23.81270  73.44995
     H 26.70726  23.81001  74.39332
     H 26.46855  20.75492  72.33939
     H 26.44825  21.40230  74.00640
     H 27.97398  21.38820  73.06968
     H 24.63622  22.94695  73.38408
     --
     0 1
     C 30.61324  20.49644  72.40614 
     C 30.47183  21.03406  71.19420 
     C 30.79531  21.59051  73.36029 
     C 30.75216  22.77934  72.62888 
     C 31.00908  21.63948  74.74382 
     N 30.43745  22.45766  71.27884 
     C 30.92671  24.02761  73.23103 
     C 31.18435  22.86580  75.34238 
     C 31.14962  24.05129  74.58925 
     H 30.33993  20.55500  70.23268
     H 31.04481  20.72825  75.32842
     H 30.86231  23.00634  70.52231
     H 30.90815  24.93911  72.64846
     H 31.36065  22.92555  76.40870
     H 31.30378  24.99933  75.08996
     H 30.61446  19.44515  72.64676

     units angstrom
     no_reorient
     symmetry c1
}

set basis 6-31g_d_p_

energy('sapt0')

10年も経つと6-31G**でもサクッと終わるのね。subaracy!

結果はこのあたりを見ればいいんですかね。ESも効いてるけどInduction(CT+mixみたいな項?)とDispersionも効いてるから、次はカチオンの有無でトリプトファンのそれぞれの原子のRESP chargeがどのくらい変化するかとか双極子どうなっているかとか確認すればいいんかな。

Special recipe for scaled SAPT0 (see Manual):
  Electrostatics sSAPT0         -13.53632077 [kcal/mol]     
  Exchange sSAPT0                10.77916633 [kcal/mol]      
  Induction sSAPT0               -5.12922298 [kcal/mol]     
  Dispersion sSAPT0              -5.69289464 [kcal/mol]     
Total sSAPT0                    -13.57927206 [kcal/mol]     
-------------------------------------------------------

Scanning the torsional potential in Psikit(RDKit+Psi4)

Considering the conformational effects of the compound is important in Structure Based Drug Design, this paper discussed about it, in terms Protein-Ligand binding using torsional scan of each ligands(PDB:2JH0,2JH5,2JH6). They calculated torsional energy, and explained the relation between inhibitory activity and torsional energy.

Torsional scanning is the task of the quantum chemistry rather than that of chemoinformatics. But I wanted to conduct quantum chemical calculation as an extention of chemoinformatics way, so I implemented it in Psikit

日本語訳

jupyterでRDKitからのB3LYP/6-31G*でのDFT計算さっくりうごくの素晴らしい。

Psi4 install error

I got a this error. Conda always irritates me ;-)

>>> import psi4
Traceback (most recent call last):
  File "/Users/kzfm/anaconda3/lib/python3.6/site-packages/psi4/__init__.py", line 55, in <module>
    from . import core
ImportError: dlopen(/Users/kzfm/anaconda3/lib/python3.6/site-packages/psi4/core.so, 2): Library not loaded: @rpath/libiomp5.dylib
  Referenced from: /Users/kzfm/anaconda3/lib/python3.6/site-packages/psi4/core.so
  Reason: image not found

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/kzfm/anaconda3/lib/python3.6/site-packages/psi4/__init__.py", line 60, in <module>
    raise ImportError("{0}".format(err))
ImportError: dlopen(/Users/kzfm/anaconda3/lib/python3.6/site-packages/psi4/core.so, 2): Library not loaded: @rpath/libiomp5.dylib
  Referenced from: /Users/kzfm/anaconda3/lib/python3.6/site-packages/psi4/core.so
  Reason: image not found

@rpathという変数がMacでは何らかの理由により正しく設定されないらしい

朝からイラツキすぎてまともに解決する気力もないので、当面放置するかも。これから買い物いかないといけないし。 それにしてもCondaよくわからんエラー多すぎるわ。RDKitもそうだったしさー。あんまよくないんじゃないの?職場では使ってないしー

forum見たけど報告されてないっぽいし、あとで気力が戻ったら投げておきます(未定)。

追記 0800くらい

エントリ書いて朝ごはん食べたら冷製になったのでなおした。

$ install_name_tool -change @rpath/libiomp5.dylib /Users/kzfm/anaconda3/pkgs/intel-openmp-2018.0.0-h8158457_8/lib/libiomp5.dylib /Users/kzfm/anaconda3/lib/python3.6/site-packages/psi4/core.so

Mishima.syk #6 やりました

参加した皆様お疲れ様でした。色々と面白い話が聞けて楽しかったです。

僕は分子設計を軽くディスってきましたw 本当にもっと頑張って欲しいですね。リードオプティマイゼーションのスタンスからFMOを使っているんだけど、もっとアーリーな段階での使い方の議論は新鮮だったので今度試してみたい。

あとは告知の仕方に問題があったので改善したい。事前に知っていれば参加したのにと何人かのヒトに言われたので。twitterやれとしか言えないが、それも難しいでしょうからメールを送るとか考えないとあかんのかな…

懇親会で興味を伺ったらディプラーニングに注目しているヒトが多いようなので次回はそれ系のネタを集めようかなと思っています。

反省会は「鈴木屋」か「やごみ」でやりましょうw

pygamessの便利なところ

懇親会で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は要らないかな。

pygamess 0.3.0

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

一番大きな変更はデフォルトでrungmsを使わないようにしたことで 直接ddikickを叩くことにしたことか。それでgamessのパスをどうしようかなぁと。

とりあえずコンストラクタで直接指定してない場合は

  1. 環境変数GAMESS_PATHを調べる
  2. 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みたいなノリで計算できるようにしたい。

python版rungms

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くらいだった。

pygamess+pybelで構造最適化計算

昨日は製薬業界の集まりがあって、他社のヒトと少し話す機会があって、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

やっぱコンストラクタに引数をわたせるようにするべきだよなぁ。

ProductName 初めてのPython 第3版
Mark Lutz
オライリージャパン / 4830円 ( 2009-02-26 )


量子化学計算したいならこれかな。

pygamessをアップデートした

Gamessのpythonラッパーであるpygamessを0.2.0にアップデートした。

主な変更点

  • pybelに対応した
  • テストをpyvowsに移行した

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はあっているなぁ。

TODO

ドキュメントをきちんと書く